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Chapter 1 
INTRODUCTION 

1.1 Motivation of Current Research 

High-speed vehicles such as the Space Shuttle Orbiter must withstand severe 
aerodynamic heating during reentry through the atmosphere. The Shuttle skin and 
substructure are constructed primarily of aluminum, which must be protected during 
reentry with a thermal protection system (TPS) from being overheated beyond the 
allowable temperature limit, so that the structural integrity is maintained for subsequent 
flights. High-temperature reusable surface insulation (HRSI), a popular choice of passive 
insulation system, typically absorbs the incoming radiative or convective heat at its 
surface and then re-radiates most of it to the atmosphere while conducting the smallest 
amount possible to the structure by virtue of its low diffusivity. 

In order to ensure a successful thermal performance of the Shuttle under a 
prescribed reentry flight profile, a preflight reentry heating thermal analysis of the Shuttle 
must be done. The surface temperature profile, the transient response of the HRSI 
interior, and the structural temperatures are all required to evaluate the functioning of the 
HRSI. Transient temperature distributions which identify the regions of high temperature 
gradients, are also required to compute the thermal loads for a structural thermal stress 
analysis. Furthermore, a nonlinear analysis is necessary to account for the temperature- 
dependent thermal properties of the HRSI as well as to model radiation losses. 

Based on the capability to handle time-dependent as well as nonlinear boundary 
conditions and thermal properties, and programmability in general purpose codes, the 
finite element method is used to discretize the governing energy equation. When the 
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structure is subjected to severe thermal loads, the discretization level must often be 
increased by adding degrees of freedom to predict accurately the temperature gradients 
and the subsequent stress response. The addition of degrees of freedom significantly 
increases the computational cost of transient nonlinear thermal analysis. Hence, it is 
desirable to employ a method which can effectively reduce the computational problem 
size while maintaining accuracy, thus enabling an efficient solution of large or complex 
thermal problems. 

The force-derivative method has shown tremendous success, in terms of the 
reduction achieved, in a variety of structural problems and also when applied to a simple 
linear transient thermal problem. This observation has motivated the present study on the 
feasibility of using the force-derivative method as a reduction technique for solving 
nonlinear transient thermal problems. 

1.2 Literature Survey 

This section summarizes the research done in the past to predict the thermal 
response with greater efficiency and accuracy in a wide range of problems. After an 
introduction to the different approximation methods available, the commonly used finite 
element time integration algorithms are discussed along with measures proposed in the 
earlier studies to overcome the demerits associated with such algorithms. Then, studies 
on the concept and usefulness of adaptive mesh generation techniques are presented. 
Finally, the evolution of the reduction methods to improve the computational efficiency 
of large-scale structural problems is outlined, followed by a review of the status of the 
reduction methods as applied to thermal problems. 

1.2. 1 Approximate Numerical Methods 

Analytical or exact methods to obtain the temperature response are often 
impossible or impractical, due to the arbitrariness or irregularity of the geometry, or other 
features of the problem. Therefore, approximate numerical methods are often employed 
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for this purpose. The thermal analysis of convectively-cooled structures by Wieting and 
Guy [1]* was based on the finite difference lumped-parameter technique such as used in 
the Martin Interactive Thermal Analyses System (MITAS) [2], The finite difference 
model of a problem gives a pointwise approximation to the governing equation of heat 
transfer. The model is formed by writing difference equations for an array of grid points, 
and hence is improved as more points are used. Bhattacharya [3] and Lick [4] used an 
improved finite difference method for time-dependent heat conduction problems. 
However, the method performed poorly when faced with irregular geometries or 
complicated boundary conditions. 

Another numerical approach is the boundary integral equation method (BEEM) 
where an exact integral formula is derived relating boundary heat flux and boundary 
temperature from a fundamental singular solution to the governing equation. The pan of 
the boundary data not already prescribed in the problem statement is obtained 
numerically from the formula. The temperature throughout the body is then generated by 
means of a Green’s type integral identity directly in terms of the boundary data. A hybrid 
method which combined the BIEM with the Laplace transform technique to solve 
transient heat conduction problems was developed by Rizzo and Shippy [5] in 1970. 
Rources and Alarcon [6] presented the formulation for a two-dimensional isotropic 
continuous solid using the BEEM with a finite difference approach in time. 

The finite element method (FEM) represented a major breakthrough in solid 
mechanics. Although the concepts of the FEM were well in use already, the method 
gained momentum in 1965 when Zienkiewicz and Cheung [7] introduced it as a method 
applicable to all field problems that can be stated in a variational form. This method 
gives a piecewise approximation to the governing equation. The finite element model 
replaces the solution region by an assemblage of discrete elements, and thereby reduces 


♦Numbers in brackets indicate References. 
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the continuum problem to one of a finite number of unknowns at points called nodes 
along element boundaries (and sometimes within the elements too). The variation of the 
field variable within the elements is expressed in terms of the nodal values of the variable 
and the assumed approximating functions called interpolating functions within each 
element. 

The finite element method was first extended to linear thermal analysis by Wilson 
and Nickell [8] and by Becker and Parr [9] to solve steady and transient heat conduction 
problems. With this introduction, the potential of the method to perform thermal analysis 
was realized as it could not only represent irregular geometry but could improve accuracy 
and sometimes could perform more efficiently for a given accuracy than finite difference 
methods. It has the added advantage that it can model both thermal and structural 
problems. The various parameters in the method that can affect the nature of the solution, 
in terms of efficiency and accuracy, have been studied extensively since the method 
evolved. ' The ease in modeling complex geometry, the capability to handle 
time-dependent as well as nonlinear boundary conditions and thermal properties, 
programmability in general purpose codes, and compatibility with a subsequent structural 
analysis have made the FEM a very useful and effective method in engineering 
applications in general. 

One of the early works in nonlinear heat transfer using finite elements was by 
Richardson and Shum [10], who included nonlinear radiation-convection heat flux 
boundary conditions in an explicit formulation. The convergence characteristics were 
improved by an alternative implicit-direct iteration scheme by Beckett and Chu [11]. 
Aguirre-Ramirez and Oden [12] applied the FEM to solve nonlinear heat conduction 
problems with temperature-dependent conductivity by the Runge-Kutta numerical 
integration scheme, while Thornton and Wieting [13] developed a procedure to handle 
several temperature-dependent parameters for simple elements based on the generalized 
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Newton-Raphson iteration technique [14]. The procedure which relies on the 
assumptions that thermal parameters are constant within an element and depend only on 
the average element temperature, was applied to convectively cooled structures [15]. 
While the early studies employed the variational principles introduced by Gurtin [16] 
(traditionally used in structural analysis) to derive the finite element equations, 
subsequent researchers used the Galerkin method of weighted residuals with interpolation 
functions as the weighting functions, commonly referred to as the conventional 
formulation. 

1.2.2 Time-Integration Algorithms 

Transient problems require the numerical solution of a set of first-order 
simultaneous ordinary differential equations. This was done by solving the incremental 
form of the governing equation [17], using an implicit time integration and modified 
Newton-Raphson iteration to establish equilibrium at every time increment. The direct 
integration techniques start from a known initial condition and then solve recursively for 
the solution at successive intervals of time based on a finite difference approximation of 
the time derivative of the temperature at an intermediate time within each time interval. 
This approximation has a significant effect on the behavior of the transient response. The 
study of the oscillation and stability characteristics of direct integration algorithms has 
received considerable attention over the years. 

In the explicit forward-difference scheme, the set of temperatures at a given time 
is expressed as an explicit function of the set of previous temperatures in the structure. It 
requires minimal computation per time step to solve uncoupled algebraic equations. 
Capacitance has to be lumped, which may have its own inaccuracies as discussed in 
reference [18]. It is only conditionally stable with severe restrictions on the time step for 
short or thin elements having high diffusivity. This makes computation costs very 
prohibitive when the response is calculated over a long duration. Implicit schemes are 
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unconditionally stable even for nonlinear problems [19] thus permitting larger step sizes. 
However, they require considerable computational effort to solve coupled equations and 
significant additional computational effort for nonlinear problems because of the need for 
iterations at each time step and the need to factor the effective coefficient matrix every 
time step. Nevertheless, implicit algorithms are more efficient for solving stiff equations 
with widely-separated eigenvalues. In the popular Crank-Nicholson implicit algorithm 
the step size is often dictated by solution accuracy, for too large a step especially in 
nonlinear applications, can introduce errors in the spatial temperature distributions. 
Myers [20] has presented a method for estimating time steps required in heat conduction 

problems. 

Adelman and Haftka [21] have identified some essential features of most transient 
heat conduction problems with respect to integration techniques, 

1. The thermal response may be divided into regions of slowly and rapidly varying 
temperatures. Steep transients accompany initial conditions or sudden changes in 

the heat load. 

2. The rapidity of variation of the transient portion of the temperature history is 
proportional to the quantity equal to the square of the length of the element 
divided by the diffusivity of the material. During such a transient, time steps 
much smaller than this quantity must be taken no matter what type of integration 
technique is used. 

3. During a period of slowly-varying temperatures, large time steps may be taken by 
implicit integration techniques but explicit algorithms still have the 
above-mentioned restriction on the time step size. 

Extensive research has been performed to improve the efficiency and accuracy of 
these algorithms. Subcycling and mixed time integration algorithms are some techniques 
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that have been investigated. Subcycling uses different time steps in different subdomains 
of a problem whereas mixed time integration uses different integrators and a single time 
step. However, application to nonlinear problems faces difficulties as the critical time 
step varies and stability characteristics are not well defined for implicit integrators. 
Orivuori [22] improved the efficiency of the Crank-Nicholson method to solve a 
nonlinear problem (nonlinearity being due to temperature-dependent material properties 
and boundary conditions), by using constant reference values for the effective coefficient 
matrices and load vector but periodically multiplying them by time-dependent functions 
to account for the nonlinearities, thus avoiding repeated factorizations. Efforts to couple 
the development of a set of various-order implicit algorithms and a strategy to 
automatically select both the largest possible time step as well as the appropriate 
algorithm throughout the solution process have resulted in the GEAR algorithms [23]. 
Reference [21] compares the GEAR algorithms against those used in SPAR [24] and 
MITAS [2]. 

The Taylor-Galerkin approach first introduced in convective transport problems 
was extended to transient nonlinear thermal-structural problems by Thornton and 
Dechaumphai [25] and applied to aerodynamically heated leading edges [26], This 
algorithm utilizes first-order Taylor series in time and Galerkin method for spatial 
discretization. Unlike conventional algorithms, this approach treats the conservation 
variable and not temperature as the unknown. The nonlinearities are conveniently 
handled through the flux components thus avoiding the need to regenerate element 
matrices for nonlinear problems. The fluxes are interpolated from nodal values, in the 
same form as the conservation variable. The resulting element matrices could be 
evaluated in closed form thus avoiding the numerical integration for complex finite 
element shapes. The merits of this linear flux approach were seen again in steady-state 
thermal -structural analysis by Pandey et al., [27]. A second-order accurate explicit 
scheme was proposed by Tamma and Namburu [28] by including higher-order time 
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derivatives in the Taylor series which are evaluated from the governing equation. An 
alternate implicit, second-order accurate approach was presented by Thornton and 
Balakrishnan [29] which uses enthalpy as the dependent variable, thereby handling 
temperature-dependent specific heat outside the element integral as well and also 
permitting larger time steps than the explicit form. 

High conductivity or very small element sizes, as may result from adaptive 
techniques, could severely restrict time step sizes from accuracy and stability 
considerations. As a means of replacing time-integration techniques and their associated 
complexities, hybrid methods have been developed which employ the Laplace transform 
technique to remove the time derivative from the governing equation and then solve the 
equation in the transform domain by the BEEM [5], finite difference [30] or FEM [31-33]. 
The temperature response is then obtained directly at the selected time of interest by 
applying an inverse transform to the solution in the transform domain. Although this 
technique was used very efficiendy for linear problems in Ref. [33], its accuracy as a 
general nonlinear solution technique is questionable. Cerro and Scotti [34] have shown 
that the linearization involved before the Laplace transform is applied, neglects the time- 
dependent behavior of the terms which define the nonlinear problem, and hence produces 
significant inaccuracies as the nonlinearity increases. 

1.2.3 Adaptive Mesh Generation Techniques 

Severe aerodynamic heating produces non-uniform temperatures and stress 
gradients over the structure and the distributions are also time-dependent. Adaptive mesh 
generation techniques are employed to capture such localized effects and thus improve 
solution accuracy for a given computational effort. The two basic approaches are : 

1. adaptive refinement/derefinement which includes 

(a) the h-method, 

(b) the p-method, 

(c) the r-method, and 
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2. remeshing. 

The h-method increases/decreases the number of degrees of freedom (DOF) by 
adding/removing elements in the original mesh in regions of interest. It is commonly 
used in production-type codes but orientation of elements cannot change to accommodate 
the varying field. The p-method maintains the geometry of the elements of the initial 
mesh but increases/decreases the DOF by increasing/decreasing the order of the 
interpolating polynomial by the use of nodeless variables. The hierarchical temperature 
interpolating functions need not be of the same order as displacement interpolating 
functions; therefore independent refinement is possible which is useful in integrated 
thermal-structural analysis [35] but has its complexities during implementation in 
computer programs. The r-method keeps the number of elements and their connectivities 
the same but relocates the nodes. This could sometimes result in distorted elements. The 
inherent drawbacks in the refinement approaches led to the use of the remeshing 
technique [36-37], wherein a new mesh is generated based on the solution obtained from 
the previous mesh. The new mesh consists of small elements in the regions with large 
changes in solution gradients and large elements in the regions where the gradient 
changes are small. For thermal problems, especially where the thermal loads move along 
the body of the structure and the magnitudes also vary with time, the mesh employed 
must adapt itself both in time and space (mesh movement) to accurately capture the 
transient temperature response. Remeshing has proved to be beneficial where steep 
gradients are involved because of high convergence rates. Mesh adaptivity based on 
relative error indicators estimated from the previous finite element solutions, avoids the 
need to know the behavior a priori. For a plate subject to nonuniform surface heating, 
when compared with a uniformly refined structured mesh, an adaptive unstructured mesh 
required much fewer nodes for a given acceptable etror or produced a much smaller error 
for a given number of degrees of freedom (nodes). 
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1.2.4 Reduced-Basis Methods 

A class of methods known as the reduced-basis methods retain the modeling 
versatility of the FEM while simultaneously reducing the number of DOF. The central 
idea is to solve the problem in a reduced subspace of the original space of discretization. 
This is done by replacing the governing equations of the structure by a reduced system of 
equations with considerably fewer unknowns. Thus, the actual solution vector is 
approximated by the Rayleigh-Ritz technique as a linear combination of a reduced set of 
linearly independent basis vectors. The approximate solution vector is then given by the 
product of a transformation matrix whose columns are the basis vectors, and a vector of 
undetermined coefficients which is obtained by solving the reduced system of equations. 

The key to an effective reduction technique is the proper choice of the basis 
vectors which may include eigenmodes, Ritz vectors, Lanczos vectors or any suitable 
combination of the above. The following guidelines aid in the choice of appropriate basis 
vectors: 

1. The vectors must be linearly independent and span the space of solutions in the 
neighborhood of the point considered on the solution path, and therefore fully 
characterize the nonlinear response in that neighborhood. 

2. Their generation should be both simple and computationally inexpensive, and 
their number can be automatically selected for any given problem. 

3. The vectors must have good approximation properties so that they provide highly 
accurate solutions on a large interval of the solution path. 

While the first property guarantees convergence of the Bubnov-Galerkin approximation, 
the other two decide the efficiency and effectiveness of the method in solving large-scale 
nonlinear thermal problems. The review of the literature available for reduced-basis 
methods is presented in separate seciiens for structural and thermal problems. 
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Structural Problems 

The use of modal methods as a reduction technique has increased by leaps and 
bounds since its introduction in 1944 by Biot and Bisplinghoff [38] to solve dynamic 
structural problems. The eigenmodes of a structure were then recognized, for the first 
time, to form a complete set of orthogonal and linearly independent vectors whose 
superposition could therefore represent the transient, linear structural response. This 
method henceforth came to be known as the mode-displacement method (MDM). Its 
attractive feature was that the reduced system of equations that resulted from this 
transformation were uncoupled, and hence could be solved individually as single degree 
of freedom systems. 

The use of modal techniques for nonlinear problems is based on the principle of 
local mode superposition. For mildly nonlinear problems, Bathe and Gracewski [39] 
successfully employed the MDM coupled with the residual force technique. Herein, a 
single set of modes (based on linear analysis) and a constant effective stiffness matrix is 
used throughout the analysis, while the nonlinearities are fully taken into account in the 
evaluation of the residual force vector. However, for highly nonlinear problems this 
could yield erroneous results as the system characteristics are continually changing, while 
an accurate solution might require too frequent basis updates which could prove to be 
expensive. Noor [40] observed that the use of the linear solution as a basis vector 
necessitates frequent additions of corrective basis vectors, where each additional vector is 
obtained by solving the full system of nonlinear equations. To avoid this, Noor suggested 
the use of a nonlinear solution and its various order path derivatives as basis vectors in 
nonlinear structural problems. To accomplish the same task, Idelsohn and Cardona 
[41-42] suggested the inclusion of derivatives of the basis eigenmodes and Ritz vectors. 
The Ritz vectors have the advantage of accounting for the spatial distribution of the load 


at the basis generation itself. 
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The efficiency of a reduction method is measured by its ability to accurately 
predict the transient response by using as few basis vectors as possible. This ability is all 
the more crucial for a modal method where the complete solution of a large eigenvalue 
problem is rather expensive and also entails huge computer storage requirements. 
Although the MDM proved to be accurate and cost effective for solving many dynamic 
structural problems, the method experienced convergence difficulties when dealing with 
problems that exhibited discontinuities in time or space, or problems with closely-spaced 
natural frequencies [50-51], The MDM often required a large number of modes to 
predict even the displacement response accurately and was ineffective in predicting the 
stresses (which is expected since the stresses are functions of the spatial derivatives of the 
displacements and the process of differentiation tends to magnify errors already existing 
in the displacement calculations). Kline [43] attempted to improve the MDM by adding a 
suitable choice of Ritz vectors to the system eigenmodes for linear dynamic problems. 
Several researchers [44-47] have worked on developing improved higher-order or faster- 
convergent modal solutions. The improved convergence of the mode-acceleration 
method (MAM) over the MDM is due to the additional term which represents the pseudo- 
static response thus including to some extent, the flexibility of the higher modes which 
are totally neglected in the MDM. 

Camarda [48-50] identified a unified approach for deriving successively 
higher-order modal methods, which are collectively called the force-derivative method 
(FDM) and offer an increasingly improved approximation of the higher modes neglected 
in the basic MDM. As the name indicates, the additional terms in the FDM involve an 
increasing order of the time derivatives of the forcing function, which are obtained by 
repeated integration by parts with respect to time, of the convolution integral form of the 
modal response. Thus the MAM, which supplements the MDM with one correction term 
that depends on the forcing function itself, is the first-order FDM. 
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The FDM has established its superiority over lower-order modal methods by the 
consistently faster convergence and greater accuracy in a wide variety of linear structural 
problems [48-50]. It must be noted that if the forcing function or its time derivatives are 
discontinuous in time, the higher-order modal methods should include appropriate jump 
conditions to avoid errors in the solution close to the time of discontinuity. This is 
because the development of higher-order methods is based on the assumption that the 
forcing function and its time derivatives are continuous. McGowan and Bostic [51] have 
demonstrated that for a multi-span beam (which has closely-spaced frequencies) subject 
to a uniformly distributed load which varies as a quintic function of time, the FDM (order 
4 and 6) not only significantly reduced the number of modes necessary to represent an 
accurate response, but also required considerably less computational time as compared to 
the lower-order modal methods and the Lanczos method. The highly desirable rapid 
convergence property of the FDM was further exhibited in the analysis of an 
unconstrained high-speed civil transport structure [52], which used an elastic flexibility 
matrix to replace the inverse of the stiffness matrix which was singular. However, a 
comparative study of the central processor unit (CPU) times in this case showed an 
increased CPU time required by the higher-order FDM as it performed more calculations 
involving the elastic flexibility matrix. 

Thermal Problems 

Transient thermal problems exhibit a wide spectrum response and the higher 
thermal modes excited by the heat supply intensities generally dominate the response. 
Bushard [53] employed the Guyan reduction technique (commonly used in structural 
dynamics) to solve transient thermal problems. The mode-superposition technique was 
introduced in thermal analysis as well by Biot [54] in 1957. However, too many DOF had 
to be retained in the reduced system for the MDM to predict an accurate temperature 
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response even for linear thermal problems. The slow convergence of the MDM can be 
traced to the omission of the important higher modes. 

Nour-Omid et al., [55] recognized that the Lanczos algorithm, which originated to 
solve the symmetric eigenvalue problem, was an efficient tool to extract 
eigenvalues/eigenmodes at both ends of the spectrum. This capability enabled the 
algorithm to generate a sequence of orthogonal vectors which served as an effective basis 
in the transient thermal analyses by Nour-Omid [56] and Coutinho et al., [57], By 
reducing the system of equations to the tridiagonal form, the solution required little 
numerical effort. The vectors were themselves inexpensive to generate. Assuming the 
spatial distribution of load does not vary with time, the steady-state solution was used as 
the starting vector for linear problems, and with few vectors generated, the coordinate 
transformation matrix contained vectors which were also global approximations of the 
higher modes [57]. Thus, Lanczos vectors served as an effective reduced-basis for 
thermal analysis. Subsequently, in an attempt to capture steep gradients in the solution, 
Cardona and Idelsohn [58] employed the increment of the nodal temperatures for the first 
time step as the starring vector to generate the orthogonal Lanczos vectors, and then 
supplemented with a constant vector for successive time steps. 

Nonlinear thermal problems were handled by Cardona and Idelsohn [58] in a 
manner similar to nonlinear structural problems by augmenting the set of basis vectors 
with derivatives of the same with respect to their own amplitude parameters, thus 
accounting for the variation of the system properties caused by changes in temperature. 
Noor et al., [59] extended the path-parameter approach to solve nonlinear steady-state 
thermal problems. The governing equation is embedded in a single or multiple-parameter 
family of equations. The path derivatives obtained by successive differentiation of the 
finite element equations of the initial discretization, are computed at zero values of the 
path parameters so that only one matrix factorization is needed. Often an augmented set 
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is used which includes a constant vector to represent a uniform temperature mode and 
reciprocal vectors along with the first few path derivatives. The cost of evaluating the 
basis vectors and generating the reduced equations can be rather high relative to the cost 
of solving the reduced nonlinear algebraic equations. This is because the expressions for 
the basis vectors grow in complexity for higher-order derivatives and their computation 
involves contractions of multidimensional arrays with the basis vectors. Besides, all the 
lower-order derivatives must be obtained before any subsequent order derivatives can be 
computed. 

On the other hand, a modified modal method reviewed in reference [40] yielded 
reasonably accurate solutions for step loaded dynamic problems by employing two sets of 
vibration modes as basis vectors, namely 

1 . Few modes of a linear eigenvalue problem based on initial conditions; 

2. Few modes of the nonlinear steady-state of the structure. 

This success prompted the application of the modified method by Shore [60-62] to obtain 
the temperature history of a model of the Space Shutde Orbiter wing subject to reentry 
heating. The nonlinearities arose from temperature-dependent properties and radiation 
from the surface. For temporally varying heat loads, provided the uniform or nonuniform 
spatial distribution of the heat loads remained constant in time, excellent results were 
achieved. This was made possible by a careful construction of the basis which included 
an adaptively generated vector based on the temperature distribution from the previous 
time interval among others. When spatially varying heat loads were considered, further 
enrichment of the basis via an analytically generated vector based on the changing 
heating distribution became necessary. 

The unified approach to develop the higher-order FDM for structural problems 
has been extended to linear thermal analysis [50]. When applied to a simple 
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one-dimensional thermal problem of a rod heated at one end, the first- and the 
second-order FDM converged faster than the MDM as they did in structural applications. 

1.3 Objective of Study 

In spite of the versatility of the FEM, an unwieldy number of degrees of freedom 
are often required to model complex geometries or to capture the temperature gradients 
arising from severe thermal loads and hence to accurately predict the stress response 
subsequently. The important step in nonlinear analyses, the solution of the system of 
algebraic equations associated with the finite element model, therefore remains very 
expensive even with improved numerical and programming techniques. This is 
especially true when analyzing large, complex structures under severe thermal effects 
since the analyses need to be carried out for long durations and involve full-size system 
matrices which result in prohibitive computer storage and analysis costs. 

Literature review indicates that reduced-basis methods have been used extensively 
to provide very efficient solutions to large-scale structural problems, and to some extent 
even for thermal problems. The modal methods use a reduced set of the lower 
eigenmodes directly as the basis vectors. The effectiveness of the FDM, measured by the 
reduction achieved in a multitude of structural problems as well as in a linear thermal 
problem, is very impressive. Nevertheless, the method has heretofore never been applied 
to nonlinear thermal analysis. The FDM has therefore been chosen as the reduction 
technique in this research effort. 

The primary objective of this study is to develop and validate a solution procedure 
that employs the FDM to obtain the transient response for nonlinear thermal problems. 
The specific objective of this study is to compare numerically the efficiency (in terms of 
the reduction in the number of equations to be solved only) of the FDM with that of the 
fundamental modal method, the MDM, in obtaining the nonlinear transient thermal 
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response (within the desired accuracy) of a realistic structure such as the lower surface of 
a segment of the Shuttle wing. 

1.4 Scope of Study 

To achieve the aforementioned objectives, first of all a new finite element 
algorithm has been developed for solving nonlinear transient thermal problems, which 
incorporates the modal methods (ranging from the MDM to the second-order FDM) and a 
fixed-point iteration scheme. The modal methods have been derived in a form that can be 
easily implemented in existing computer programs. This fact coupled with the desire to 
take advantage of existing advanced finite element software, has led to the 
implementation of the new algorithm in the Computational MEchanics Testbed 
(COMET) system [63]. 

The analytical solution has been used to solve the reduced system of uncoupled 
modal equations, thus eliminating the need to employ a finite difference approximation in 
time to solve the full-system of finite element equations. The solution, which is in the 
form of a convolution integral, is obtained by stepping in time though, in order to 
facilitate a piecewise linear approximation to all nonlinear quantities involved, thus 
minimizing the error that could result from a totally linearized approach. In this 
approach, the obvious restriction on the time step size is imposed by the degree of 
nonlinearity of the problem; unlike in conventional time-integration algorithms, where 
the time step size is directly dictated by the stability criteria. A study of the effect of time 
step size and frequency of eigensolution updates on the solution accuracy has been 
performed, to a limited extent, for a one-dimensional nonlinear problem. 

A basic understanding of the role played by the correction terms of the higher- 
order modal methods in predicting the response accurately with fewer modes has been 
provided. Guidelines are presented to show how the relative rate of convergence of the 
correction terms as more modes are included in the solution, can be used to make a priori 
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estimate of the number of modes required. Attempts have also been made to examine the 
effect of the finite element mesh on the efficiency of the reduction process. Finally, to 
assess the applicability and the effectiveness of the FDM as a reduction technique, the 
method has been applied to perform the thermal analysis of the lower surface of a 
segment of the Shuttle wing. 

A statement of the general transient heat conduction problem and the finite 
element formulation, both linear and nonlinear, are presented in chapter 2. For linear 
problems, the modal methods up to the second-order FDM, the modal coordinates for the 
specific case of a linearly time-varying load, and guidelines for a priori estimate of the 
required number of modes are all derived in chapter 3. A simplified Newton-Raphson 
iteration scheme, and the modified modal methods to obtain the transient response for 
nonlinear problems are described in detail in chapter 4. The results of linear example 
problems are discussed in chapter 5. Nonlinear thermal problems in one and two 
dimensions have been solved, and the results are discussed in chapter 6. A summary of 
the conclusions and suggestions for future research appear in chapter 7. 



Chapter 2 

GOVERNING EQUATION FOR TRANSIENT HEAT CONDUCTION 

In this chapter, a finite element formulation is presented for the computation of 
transient temperature distribution in solids with general surface heat transfer. It begins 
with a statement of the law of conservation of energy which considers the work done by 
the stresses, thermal energy transported across surfaces by conduction, thermal and 
mechanical energies stored within the material, and kinetic energy due to deformation. 


The energy equation for a continuum in solid mechanics is 


acn 

9xi 


d g U 

at 



Qint 


( 2 . 1 ) 


where the subscripts i and j are dummy summation indices ranging from 1 to 3, qi are the 
heat flux components in the coordinate directions xi, Ojj are the stress components, £jj are 
the strain components, u is the unit internal energy, p is the density, and Qint is the 
volumetric rate of internal heat generation. The equation indicates a thermomechanical 
coupling which is the conversion of mechanical to thermal energy. Extensive studies on 
coupled thermoelasticity have shown that this coupling effect can often be neglected in 
the analyses of flight structures, because the thermal energy converted from mechanical 
energy is insignificant compared to the external energy resulting from intensive 
aerodynamic heating [64]. Based on this assumption, the solid is considered 


undeformable and the internal energy becomes a function of temperature alone. To be 
specific, the term cry may be dropped from Eq. (2.1) and can be expressed as 
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2 1 Problem Statement 


With the above-mentioned simplified 011 ’ the energy equation in Cartesian 
coordinates for a general three-dimensional anisotropic solid of volume V bounded by a 
surface S, Fig. 2.1, is 


§9i 

3xj 


3 (cT) 


+ Qint = P 


( 2 . 2 ) 


where p is the density, c is the specific heat, and the heat flux components are given by 
Fourier’s law as 


3i = 



(2.3) 


The material properties p, c and kjj may be temperature-dependent, where ky are the 
components of the symmetric conductivity tensor, K Substituting Eq. (2.3) in Eq. (2.2) 
we obtain the governing heat conduction equation which is solved subject to an initial 
condition 

T (xi, X 2 , X 3 , 0) = To (xi, X 2 , X 3 ) (2.4) 


and general boundary conditions which include prescribed boundary temperatures T s on 
surface S 1 , specified surface heating on S 2 , convective heat exchange on S 3 , and incident 
and/or emitted radiation on S 4 


T s = Ti (Xj, x 2 , x 3 , t) 

on Si 

(2.5a) 

qi ni = - qs 

on S2 

(2.5b) 

qi ni = h ( T s - T$) 

on S 3 

(2.5c) 

qi nj » oeT s 4 - aqr 

on S4 

(2.5d) 


The prescribed temperature, Ti, may be a function of position and time, the specified 
heating rate, q s , could be time-dependent, the convective heat transfer coefficient, h, may 
depend on the environment temperature, T e , and/or time, T s is the unknown surface 
temperature, o is the Stefan-Boltzmann constant, £ is the surface emissivity that may be a 
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function of T s , a is the surface absorptivity, q r is the rate of incident radiation per unit 
area, and ni are the components of an outward unit normal vector. 

2.2 Finite Element Formulation 

The solution domain is discretized into elements with r nodes each. The 
conventional approach is to express the temperature and temperature gradients within an 
element in terms of the interpolation functions and their gradients. The element 
temperature is defined as 

T e (xi, X2, X3, t) » N (xi, X2, X3) T (t) ( 2 . 6 ) 


where N is the row vector of temperature interpolation functions and T (t) is the column 
vector of nodal temperatures. The components of the temperature-gradient inteipolation 
matrix are given by 


5N 


By (xi, X2, X3) = i=l, 2 , 3 ; j*l, 2 , ...r 


(2.7) 


Starting with the energy equation, applying the method of weighted residuals, integrating 
using Gauss Divergence Theorem, and finally introducing the boundary conditions results 
in a set of nonlinear finite element equations given in matrix form as 

C(T) t + (Kc(T) + K h (T, t) + Kr(D ) T(t) 

= Rq (T, t) + Rq (T, t) + R h (T, t) + R r (T, t) (2.8) 


The details of the manipulation required to obtain Eq. (2.8) can be found in Ref. [18], 


The element capacitance matrix C, and the coefficient matrices K c and Kh related to 
conduction and convection respectively, are defined in Eqs. 2.9a-2.9c below, while the 
radiation matrix, K r , implicitly given by Eq. (2.9d) is explicidy given in Eq. (2.9e): 



pc N T N dV 
B T KBdV 


(2.9a) 


(2.9b) 
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K h = 1 h N T N dS 

Js 3 

(2.9c) 

K r T = I GET 4 N T dS 

Js 4 

(2.9d) 

K r =| aeT 3 N T NdS 

JS 4 

(2.9e) 

The heat load vectors due to internal heat generation, specified surface heating, surface 
convection and incident surface radiation, Rq, R q , Rh, and R r respectively, are defined as: 

Rq * ^ Qint N T dV (2.10a) 

Rq = 1 q s N T dS 

Js 2 

(2.10b) 

R h = 1 h T e N T dS 

JS 3 

(2.10c) 

R r = j aq r N T dS 

JSa 

(2.10d) 


For the sake of brevity, Eq. (2.8) is hereafter written as 

C(T) t + K(T,t) T(t) = R(T,t) (2.11) 

where K (T,t) is the system conductance matrix which contains the contributions from 
Eqs. (2.9b) - (2.9d) and R (T.t) is the system combined load vector made up of the 
vectors defined in Eqs. (2.10a) - (2.10d). Equation (2.11) is a general nonlinear 
formulation of finite element equations for transient heat conduction in an anisotropic 
medium, and the solution requires an iterative scheme combined with a suitable time- 
integration scheme. For linear thermal problems K is independent of temperature. If K 
is further simplified to be time-independent, Eq. (2.1 1) becomes 

Cf + KT(t) = R(t) (2.12) 

The solution of Eq. (2.12) requires a time-integration scheme alone. 


Chapter 3 

THE FORCE-DERIVATIVE METHOD FOR 
LINEAR TRANSIENT THERMAL PROBLEMS 

3.1 Unified Derivation of the Modal Methods 

For the purpose of completeness, the mode -displacement method used to obtain 
the transient response for linear thermal problems is first presented, followed by the 
unified derivation of the higher-order force-derivative method assuming the forcing 
function possesses continuous derivatives. 

3.1.1 The Mode-Displacement Method (MDM) 

The governing finite element equation for a linear transient thermal problem, 
Eq. (2.12), is reproduced below 

Ct + KT(t) = R(t) (3.1) 

with initial condition 

T(0) = To (3-2) 

where in general, the applied load vector R, may be time-dependent. As mentioned at the 
end of chapter 2, although the system matrix K can be time-dependent for linear 
problems, it is assumed constant in this study. The solution to the homogeneous form of 
the linear conduction equation Eq. (3.1) is given as 

T(t) * e'*** 0r; r - 1, 2. . . n (3 3) 

where <J> r is a modal vector of unknown amplitude, X r is the associated decay constant 
(analogous to the natural frequency in structural dynamics), r is a summation index and n 
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is the number of unconstrained degrees of freedom. Use of Eq. (3.3) in Eq. (3.1) results 
in the following constant eigenvalue problem (EVP): 

K <» r - X r C <fr r = 0 ; r = 1, 2. . . n (3.4) 

where n is the total number of degrees of freedom. The eigenvectors are normalized such 
that 

$ p C <j) q = 5pq; p,q = 1, 2. . . n (3.5) 

and 

<()p K <j)q = 5pq ^pq» p>q = 1 , 2 . . . n (3.6) 

where 8pq is the Kronecker delta. In matrix form (Eqs. (3.5) and (3.6) become, 

d> T C <D = I (3.7) 

and 

<D T KO = A (3.8) 

where I is the identity matrix and A is a diagonal matrix with entries that are the 
eigenvalues. The homogeneous solution to Eq. (3.1) in the form of a modal summation is 
given by 

T(t) = X 4> r zr(t) (3.9) 

i*i 

where the solution is expressed as a linear combination of all the eigenvectors of the 

system weighted by the unknown modal coordinates. 

Substitution of Eq. (3.9) in Eq. (3.1), premultiplication by <|>r . and the use of the 

orthogonality of the modes, Eqs. (3.5) and (3.6), result in the following uncoupled modal 


equations: 


T 

z r (t) + Xf Zp(t) * 4>r E(t), r = 1 , 2. . . n 


( 3 . 10 ) 
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Pre-multiplying Eq. (3.9) by <t> r T C and setting t = 0, the initial condition becomes 

ZrO = <t>r CTo (3.1 1) 

The analytical solution to the nonhomogeneous problem represented by Eqs. (3.10) and 
(3.11) is given as follows: 

z r (t) -^c-M + f e^frRWdx (3.12) 

JO 

T 

where the integral is called the convolution of the functions e***© and <)) r R(t). 

Approximating the solution to Eq. (3.1) by using a truncated set of modes, one has. 

P 

T(t) 2 X <|>r z r (t); p*n (3-13) 

r=*l 

where z r (t) is given by Eq. (3.12) and <t>r» and X r are obtained from the solution of 
Eq. (3.4). Rewriting Eq. (3.13) in matrix form results in 

T(t)sSz(t) (3.14) 

where A denotes a reduced set of modes. The approximate solution given by Eq. (3.13) or 
(3.14) is commonly referred to as the mode-displacement method (MDM) in structural 

dynamics. 

3.1.2 The Mode- Acceleration Method (MAM) 

Equation (3.10) may be rearranged as follows: 

Zr (t) = -L<|fR(t)-- L z r (t) (3.15) 

Xf Xf 

T(t) - X R(<> - "L 

£ »1 *• 


So Eq. (3.9) becomes 


(3.16) 
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Using the orthogonality of the modes and taking the inverse of the square matrix on either 
side of Eq. (3.8) yields 

d>W T = ± (3-17) 

A 

T 

Pre- and post-multiplying both sides of the above equation by 0 and O , and noting that 
d> 0 1 = O *0 = I , leads to the following identity: 

K' 1 =0-L0 T (3.18) 

A 

or 

4>r T ( 3 ‘ 19 ) 

r=l ** 

Now Eq. (3.16) may be written as follows: 

T(t) = K' 1 R(t) • X 4>r z/O (320) 

i»l ** 

8f (x*cp) 

The Leibnitz’s rule for differentiation of an integral states that if f (x,<p) and — are 


continuous functions of X and <p, then 


jL 

dtp 


J 'V(tp) 
u(<p) 


f(x.<p) <*x 


p ,<<p) 3f(2t2l (ix . iLf[|i(9),<p] + (3.21) 

L«p> ^ d "> d,,> 


provided p. and v have a continuous first-order derivative with respect to <p. If the forcing 
function is C° continuous, then the Leibnitz’s rule can be applied to Eq. (3.12) to replace 
the term z r (t) in Eq. (3.20) as follows: 

,J e W ' T) <|J R(t) dt (3.22) 

Jo 


zr(t) = 4*R(t)-^z I oe Xft 


Substituting Eq. (3.22) in Eq. (3.20) and rearranging gives: 

T(t) = X 4>r Zr(t) + (k" 1 - X ^ft T ] R W 


r=l 


(3.23) 
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Now introduce an approximation to the solution by using only a subset of the modes, 
so that 


p 

T(t) = X <J>rZr(t) 


r=l 



(3.24) 


or, in matrix notation 

T(t) = <j>Z(t) + |k _1 -8^3 JrW (3-25) 

Equations (3.24) and (3.25) are alternate forms of the MAM, where the second term 
represents the first-order correction to the approximate solution given by the MDM, 
Eqs. (3.13) and (3.14). 


3.1.3 The Force-Derivative Method 

Assuming sufficient regularity, Eq. (3.15) may be differentiated once with respect 
to time to obtain 

z r (t) = -L<ljR(t)-- L z r (t) (3-26) 


Substituting the above in Eq. (3.15) results in the following equation: 


T X 

Zf (t) = -L<t>r R(t) - -Lfc R(0 + -Vzr(t) 
^ £ £ 


(3.27) 


So Eq. (3.9) now becomes 

T(t) = y ** + ^ W (3 ’ 28) 

i*l ** -1 V M 

Again, based on the orthogonality of the modes, the following identity can be established: 

K^CK' 1 *^ J“<D T (3.29) 

A 2 


or 


K 1 C K 1 


i-i V 


(3.30) 
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so that 


T(t) = K' 1 R(t) - K^CK R(t) +X <t>r Zr(0 

X r 


(3.31) 


If the forcing function is C 1 continuous, employ the Leibnitz’s rule for differentiation, 
Eq. (3.21), to evaluate z r (t) given in Eq. (3.31) results in 


z r (t) = -X r <jJ R(t) + R(t) + X T zrf) e'V + X? I e ^ l ' x) <(J R(t) 

JO 


dt (3.32) 


Using Eq. (3.32) in Eq. (3.31) and rearranging results in 


T(t) = £ <t»r zr(0 + [ K X f <l>r T ) R (0 + 1 - K 4 C K' + X d> r <Dr‘ R (t) (3.33) 


r=l V r=l / \ r=l X r / 

Again, approximating the response by using a subset of the modes for the summation 
terms in Eq. (3.33) one obtains: 

T(t) = X 4>r Zr(t) + f K 4 - X r 1 " <l>r T ) R (0 

ml V n«l ** ) 


+ |-K CK +X<t> r -L «|>r j R(t) 


or, in matrix form 


TtOsOZfO+jK* 1 jR(t) 


-K^CK^+S-l-S fc(t) (3.35) 

A I 

Equation (3.34) or (3.35) is the response by the second-order FDM, where the third term 
represents the second-order correction to the approximate solution given by the MAM, 
Eqs. (3.24) and (3.25). 




In general, for a O’ continuous forcing function, the highest (j+l)th order FDM for 
approximating the response is given as follows: 



^ ^ -m 

+ O (-A) 



R 


(m-l) 


(t) 


(3.36) 


where the superscript in parentheses of R(t) denotes the order of differentiation with 
respect to time. 

The major advantage in using the analytical solution, Eq. (3.12), of the uncoupled 
modal equations is that when the load vector is available as a function of time such that 
exact integration is possible, the solution at any time can be obtained directly with time 
t = 0 as the initial condition. Furthermore, the only source of error is the truncation of 
modes which can be alleviated by the proper selection of modes. However, when the 
heat load vector is a function of time that cannot be integrated exactly, numerical time 
integration is required and the time step size for a desired accuracy will depend on the 
order of the numerical integration method used. The following section describes how to 
obtain the modal coordinates for a specific type of load. 


3.2 Modal Coordinates for a Linearly Time-Varying Load 
If the applied load varies linearly with time and has a constant slope over the 
entire time domain as shown in Fig. 3.1, then the modal coordinates given by Eq. (3.12) 
can be obtained by exact integration. Rewriting Eq. 3.12 for this case one obtains: 

^(t)-z I oc‘ X,t + f e* XT(t ' t) <Dr T R(t)dt (3.37) 

JO 

where 

R (t) = R (0) + R (0) r, 0 £ t £ t (3-38) 

It is clear from Fig. (3.1) that R (0) is the inital load vector and R (0) represents 
the constant slope. Substituting Eq. (3.38) in Eq. (3.37) and expanding results in: 
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Fig. 3.1 Forcing function that varies linearly with time. 
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z^(t) = Zioe’^ 1 + 4>r R (0) l e^ tx ^dx 

JO 

+ <t> r T R(0) [ te'^^ dx (3.39) 

JO 

Using integration by parts and regrouping the terms, one has 

z r (t) = z r0 e' Xrt + ^- [r (0) ( 1 - e~ V 1 + R (0) (t • 1- s^ -) | (3-40) 

Now consider a more general case where the applied load vector can be treated as 
a series of piecewise linear functions of time with varying slopes over the entire time 
domain. Then Eq. (3.40) can be applied individually to each time interval, dt, using the 
modal coordinates obtained at the end of the previous time interval as the initial 
condition. Accordingly, Eq. (3.37) may be modified as follows: 

J fd 

e *X,(dt-x)^J R(t)dt (3.41) 

0 

where R (x) = Rn-l + Rn-ix» O^x^dt (3.42) 


Similarly, it can be shown that eq. (3.40) now becomes 


Zn.r = 2n-l,re* Xf<Jt +<t> 





(3.43) 


3.3 A Priori Estimate of the Required Number of Modes 
To minimize the computational costs associated with the complete solution of the 
EVP, some guidelines are presented below on selecting the number of modes prior to the 
actual solution process. These guidelines are limited to linear transient thermal problems. 

The error introduced in solving Eq. (3.1) by the MDM, Eq. (3.13), can be 
evaluated based on the level of approximation of the heat load vector represented by that 
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subset of modes [65]. The participation factor of the rth mode, rj r , at time, t, is defined as 
follows: 


Tlr(t ) = <t) r T R(t) (3.44) 

which gives a measure of the significance of that mode in representing the total load and 
hence the response. Pre-multiplying Eq. (3.44) by C<|> r , and using the orthogonal property 
of the modes with respect to C, the approximate load vector is given by 

P 

R(t)s£ C<Mr(t); P< n (3-45) 

r=l 


Once a desired level of approximation of the load vector has been attained, the generation 
of modes may be terminated. The error, ei, in the load vector represented by the 
truncated set of modes is given by 


/ R(t) T (r (t) - R(t)) 
V R(t) T R(t) 


(3.46) 


If a direct relationship between the errors in the solution and the load vectors can be 
established then Eq. (3.46) may be used to obtain a priori estimate of the number of 
modes required by the MDM to achieve a desired degree of solution accuracy. 

To determine the truncation of modes for the higher-order methods, which has not 
been studied in the past, the corrections offered by these methods to the MDM is 

considered. The full-system solution may be split into two parts given as follows: 

P n 

T(t) = £ 4> r z/0 + £ <t> r Zr(t) (3.47) 

i*l r*p+l 

where the first term is the approximate solution given by the MDM, and the second term 
represents the contribution of the higher modes neglected by the MDM. As before, using 
the orthogonality conditions and assuming the transient load vector has continuous 
derivatives up to order j, it can be shown that 
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" n n t n 

X z r (0 = X r^- <\>J R(t) - X 0r 0r R(t) + X <i>r 4>r R(0 

r=p+l r =p+l -*-r r=p-+-l X r r=p-t-l x/ 

-••• + X <-l) j 0 r — 1 — R^(0 + X (- 1 ) j+1< t>r — L-r z^ J+1 ^(t) (3.48) 
^ (X r ) j+1 "*« M 

For the specific case where the load vector has non-zero derivatives up to first-order only, 
that is, j = 1 in Eq. (3.48), neglecting the last term in Eq. (3.48), one has: 

X <i>rz r (t)= X 4>rf <Dr T R(0- X 4>r-^$?R(t) (3-49) 

r=p+l r=p+l r=p+l X r 

It is easily recognized that the first term on the right-hand side of Eq. (3.49) is the 
correction offered by the MAM, Eq. (3.24), and will hereafter be referred to as CMAM 
for brevity. That is 

CMAM* X - (3.50) 

rspfl 

or 

1 P T 

CMAM = K' RW-Xfcr 1 - <l>r R(t) (3.51) 

r»l ** 

Thus, by adding the contribution of the higher inodes as well to the pseudo steady- state 
response, the MAM yields a better approximation of the total transient response than 
the MDM. 

In a given spectrum, the magnitude of the eigenvalues successively increases. 
Examination of Eq. (3.50) indicates that at a given time, the significance of CMAM 
should decrease when higher modes are included in the MDM, as the eigenvalue appears 
in the denominator. In other words, the lower modes play a more important role in 
approximating the pseudo steady-state response. Based on this fact, the number of modes 
required to realize the maximum benefits of the MAM can be determined when the 
pseudo steady-state response has been approximated within a desired accuracy. 
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Alternatively, similar to Eq. (3.46), a measure of the error, e 2 , in the representation of 
K' 1 R (t) by the lower modes can be obtained as follows: 

(3.52) 

The second term on the right-hand side of Eq. (3.49) can be identified as the 
correction offered by the second-order FDM, Eq. (3.34), to the MAM and will be referred 
to as CFDM for brevity, that is, 

CFDM = - X 4>r-^ <t> J R (t) (3.53) 

r=p+l X r 

or 

. P _ 

CFDM = - K CK' + (0 (3 54) 

r=l Xr 

By including the contribution of the higher modes in approximating the first-order 
transient response, - K' 1 CK* 1 R (t), derived from the pseudo steady-state solution, the 
FDM enhances the MAM. At a given time, the magnitude of CFDM declines at a faster 
rate than that of CMAM since the square of the eigenvalue appears in the denominator. 
This makes the selection of modes more and more distinct as the order of the method 
increases. Similar to the MAM, the number of modes necessary for the convergence of 
the FDM is determined when - K' 1 CK* 1 R (t) is approximated by the lower modes 
within a specified tolerance, or the error norm, e 3 , given below becomes acceptable: 

(3.55) 

Although higher-order corrections involve the eigenvalues raised to successively 
higher negative exponents, they need not necessarily be negligible compared to the 
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lower-order corrections at all times, because they are dependent on the transient nature of 
the load as well. With just one mode included in the solution, a comparison of the norm 
of CMAM and CFDM, Eqs. (3.51) and (3.54) respectively, indicates their relative 
importance at that particular time. This rate of convergence with respect to modes 
depends on the nature of the variation of the eigenvalues within the spectrum. For a 
linear thermal problem subject to a specified transient load, it is thus possible to select the 
higher-order method which can achieve the maximum reduction at any given time. 



Chapter 4 

METHOD OF SOLUTION FOR 
NONLINEAR TRANSIENT THERMAL PROBLEMS 

4.1 Linearization of the System of Equations 
After a brief introduction to the Newton-Raphson scheme, its application to the 
nonlinear, transient heat transfer equation is described and the simplifying assumptions 
made in this study are highlighted. Although the final linearized form can be obtained by 
a direct application of the fixed-point iteration scheme, the details are presented here to 
aid in the future use of a rigorous form of the Newton-Raphson scheme. 


4.1.1 Derivation of the Newton-Raphson Method 

Consider a typical nonlinear set of equations of the form 

F (T) = 0 ( 4 - ! ) 

If T is the exact solution vector, then Eq. (4.1) is satisfied identically. However, often 
one cannot compute the exact solution, but can obtain an approximate solution so that the 
unbalance in a typical equation is smaller than a specified tolerance. This is achieved by 
the Newton-Raphson scheme which is derived below by an intuitive approach based upon 

the Taylor polynomial [66]. 

Suppose that F is twice continuously differentiable in the interval [a,b]. Let 
T e [a,b] be an approximation to T such that F'(T) is non-singular and the difference 
between TandT is small. Considering the second-degree polynomials for F(T), 
expanded about T one has: 
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— T 

F(T) = F(T) + F'(T) (T-T) + F"(^(T))(T-T) (4-2) 

where ^(T) lies berween T and T. But F(T) = 0, and since the difference between T and T 
is assumed to be small, square of the difference is even smaller and can be assumed 
negligible. Then, 

0-F(T) + F'(T)(T-T) (4.3) 

which gives 

T-T-F'OVVcr) (4.4) 

which is a better approximation to T than T. This sets the stage for the Newton-Raphson 
method, which involves generating the sequence of iterates T defined by 

T‘ = T*" 1 - FTT 1 ' 1 ) W*' 1 ), i > 1 (4.5) 


which may be rewritten as 


j 1 * 1 at' = -f i1 

(4.6) 

where 



(4.7) 

is called the Jacobian matrix and 


_i _i-l 
AT = T -T 

(4.8) 


It is clear that Newton’s method cannot be continued if the Jacobian is singular 
for some i. The method is most effective when J is bounded away from zero near the 
fixed point . Although Newton’s method will sometimes converge even with a very poor 
initial approximation, in many cases it is imperative that a good initial approximation be 
chosen. Also, Newton’s method will converge quadratically under suitable conditions 
onF. In fact, Newton’s method can be derived as a special case of the fixed-point 
iteration scheme which exhibits linear convergence in general. 
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In practice, the method is generally performed in a two-step manner. First a 
vector AT* is found which satisfies Eq. (4.6). Once this is accomplished, the new 
approximation T is obtained from Eq. (4.8). A difficulty in this method arises from the 
necessity to invert the Jacobian matrix during each iteration which involves reforming 
and factorization at each iteration, making the method computationally expensive for 
large sets of equations. In a modified form of the method, the Jacobian is formed and 
factored only once and is held constant throughout the balance of the iteration process 
[18]. More iterations are required with the modified method, but usually net 
computational costs are reduced. 


4. 1 .2 Application of the Newton-Raphson Method with Simplifications 

The discrete nonlinear transient heat conduction equation at any time t, Eq. (2.1 1) 
is reproduced below: 

C(T)f +K(T,t)T(t) = R(T,t) (4.9) 


The residual or unbalanced load in the nonlinear system of equations at the (i-l)th 
iteration is given by, 


F(t l ’\ T 1 ' 1 ) = CCT 1 ' 1 ) f ‘ * + KCr 1_1 .t) T w (t) - R(T"\t) 


ii’1. 


i-1 ~i-l. 


,i-l 


(4.10) 


Simplifying the notation, a typical pth equation of this system is given by, 

Fp 1 = Cpri tin 1 + Kpm(t) T^(t) - R‘‘ l (t); p = 1, 2 . • . n (4.1 1) 

where m is a dummy summation index, and n is the total number of equations. The 
sequence of successive iterates generated by Newton’s method for the temperature and its 


time derivative increments is then given by 


Jlps ATs + J2ps AT s = 


.i-l 


(4.12) 


where Ji and J 2 are the Jacobians associated with temperature and its time derivative, 
respectively, and are defined according to Eq. (4.7) as follows: 
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(fj^T M “ dj2PS U)t“ 


The new values of the unknowns at the end of the current iteration will then be 


„i .J-l 


i . i-1 


T s = Ts + AT S and T s = T s + AT S 


(4.14) 


While the right hand side of Eq. (4.12) is given directly in Eq. (4.10), the product terms 
involving the Jacobians on the left hand side of Eq. (4.12) are derived from Eqs. (4.13) 
and (4.10) using the chain rule of differentiation as follows 


Jips 1 ATs = Cpm (H*) Afj = c£ 8^ ATs = c£ ATs 


(4.15) 


T AT* ^Cpi 

J 2ps AT S = 


Trf) "Vm AT S ' ♦ C (|!^) "at; 




Equation (4.12) now becomes 

Cps 1 At j + [ ACps + Kps + AK ps - AR ps ] M At| = -Fp" 1 (4. 17) 

The incremental temperature-dependent system matrices and load vector are given by 



(4.18) 


(4.19) 


(4.20) 


Dropping the indicial notation, Eq. (4.17) may be rewritten as follows: 

C M At 1 + [AC + K + AK - AR f 1 AT* = - F*" 1 (4.21) 

It is seen that the Jacobian associated with temperature involves the increments in 
the system nonlinear capacitance, conductance matrices and load vector, namely AC, AK 
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and AR. These increments reflect the change in material properties with temperature, and 
generally have a minor effect on the accuracy of the Jacobian J 2 [24], The matrix AR 
includes the effect of radiation exchange as well. The evaluation of these incremental 
matrices is cumbersome and prohibitively expensive. Even though the system matrices C 
and K are symmetric their increments are not, and AR is unsymmetric as well. 
Consequently, the iterative solution requires an unsymmetric solver. Moreover, 
evaluation of the true Jacobian requires additional knowledge about the nature of 
nonlinearities and identification of element types. Therefore, for ease of implementation, 
the Jacobian is approximated by the conductance matrix alone neglecting the incremental 
quantities, which is exact if material properties are not temperature-dependent and there 
is no radiation exchange. This makes computation of the Jacobian, as well as the solution 
of the linearized system of equations easier as matrix symmetry is maintained, albeit 
convergence could be slower. Introducing this modification in Eq. (4.21), we now have 

C 1 ' 1 At 1 + K 1 " 1 AT* * -F*' 1 (4.22) 

where 

K = Kc + Kh+4K r (4.23) 

Substituting for AT and AT* from Eq.(4.14), using Eq. (4.10) and simplifying, one 
obtains the linearized equations given as follows: 

C i * 1 t‘ + K l * l (t)T*(t) = R l ‘V) (4.24) 

where the temperatures at the current iteration, and not their increments, are the 
unknowns. 

4.2 The Force-Derivative Method for the Transient Response 
4.2.1 One-Step Approach 

For constant or linearly varying transient loads, the approach used to obtain the 
transient response for linear problems in chapter 3 can be readily extended to solve 
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Eq. (4.24). Thereby, the response is computed at the time of interest in one step, with 
time t = 0 as the initial condition. The initial eigenvalue problem of Eq. (4.24) is given 
by, 


KoOo-AoCo<I>o = 0 (4.25) 

where 

Ko = K(T 0 ) and Co = C(T 0 ) (4.26) 

However, due to the nonlinear nature of the system matrices, the eigensolution is varying 
with time as well. The eigenvalue problem at time t is 


i-1 i-1 i-1 _i-l _ i-1 A 

K d> -A C <D =0 


(4.27) 


which needs to be updated during the iteration process. 


Assume a solution to Eq. (4.24) in the form of a modal summation 


T*(t)*]£ (t^zkt) 
M 


(4.28) 


which yields the uncoupled modal equations given as follows: 

T . 

if(t) + X r i * 1 z r i (t) = <|>;* 1 R ll (t) 


(4.29) 


with initial condition 


z r (0) * Zj<> =4>iO CqTo 


(4.30) 


The solution to Eq. (4.29) analogous to Eq. (3.12) is given by, 

't 

Zr(t) ■ ZK>e‘^ f0t + I e*^ il(l * T) <|>r 1 R(T l l (t),x) dt 


(4.31) 


It is important to note that for nonlinear problems, the load vector and the system 
matrices, and hence the eigensolution, are functions of temperature and are therefore 
continually changing as the temperature distribution evolves with time. Careful 
inspection of Eq. (4.31) reveals that a major approximation is embedded in that equation. 
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Since the solution is not computed at intermediate times, the nonlinear quantities are 
evaluated directly as a function of T (t), the current estimate of the temperature at the 
desired time of response. This essentially ignores the variation of the integrand with time 
due to nonlinearity, and includes only the explicit time-dependence when evaluating the 
integral in Eq. (4.31). Such an approximation could result in large errors for highly 
nonlinear problems. It is therefore imperative to introduce a time-stepping scheme so that 
the transient variation of the nonlinear quantities are well represented when obtaining the 
modal coordinates. If the transient variation of the load is of a complex nature, the 
one-step approach is not applicable in the first place, similar to linear problems. 

4.2.2 Multi-Step Approach 

Since Eq. (4.24) holds true at any instant of time, the temperature Tj at the current 
iteration at time t n must satisfy the following equation 

C„ 1 1„ +Kn' 1 Ti = Ri* 1 (4.32) 

where, the subscript n denotes the current computation time, and the superscript i denotes 
the current iteration number, and the following simplifying notations are used, namely 

kJ* = K (Tn* 1 , t„) (4.33) 

C^-Ccrf* 1 ) (4.34) 

Rn 1 » R (Tn* 1 , tn) (4.35) 

The solution is marched out in time from the initial temperature at time t = 0, and the time 
step is defined as follows: 

dt = t„ - t n -i (4.36) 

The initial condition to start the iteration process at each time step is given by, 

Tj.Tn.iJ i= 1 


(4.37) 



44 


Although the temperature-dependent thermal properties and load vector are allowed to 
vary during the iterations at each time step, the eigensolution, which is determined from 
the converged temperature at the end of the previous time step, is held constant until 
convergence at the current time step. To accomodate this piecewise linear approximation 
of the eigensolution, the left-hand side of Eq. (4.32) is modified and the change is moved 
over to the right-hand side to form a corrective heat load vector Qnl^ 1 - Accordingly, we 

now have 

C n -l ti + K„.i Tn = Rn 1 + Qn^' 1 (4.38) 

where 

= - ( [Kn* 1 - K„-i] T i' 1 + [C„ 1 - C n -i] T„ S (4.39) 

In a more general approach, Eqs. (4.38) and (4.39) may be rewritten as follows: 

C n . k ti + K n .k T‘ * R [ n l + QnC 1 (4.40) 

and 

QmZ 1 = - ( [Kn* 1 - Kn-k] Tn 1 + [CA" 1 - C n -k] Tn * ) (4.41) 

where k is an integer greater than zero. Here the eigensolution is held constant over a 
period of several time steps, say m, depending on the severity of the nonlinearity; 

k indicates the number of time steps since the previous update. The applied load vector 
may be split into a linear component, Rl„, (which may be time-dependent) and a 

nonlinear component Rj^' 1 . The entire right-hand side of Eq. (4.40) may be considered 

as a generalized heat load vector which is defined as follows: 

Q (Tn \ tn) * Qn 1 ■ Rl„ + R NlA' 1 + QniA * 1 ( 4 42) 

Another major variation between the linear and nonlinear solutions is in the definition of 
the right-hand side load vector Eq. (4.42). For linear problems, it comprises the linear 
applied load vector alone, while it is rather complicated for nonlinear applications where 
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a nonlinear applied load vector may exist and additional corrective nonlinear terms are 
involved as well. The eigenvalue problem associated with Eq. (4.40) is given by, 

Kn-k ^n-k ' A n -k C n -k ^n-k = 0 (4.43) 


The modal equations now become 

. i * i a T -.i-l 
Zn j + ^n-k,r z nj = 9 n-k,r Qn 


(4.44) 


with initial condition 


z n-lj = 9n-k/ Cn-k T n -l 


(4.45) 


The solution to Eq. (4.44) is given as follows: 


J fdt 

e *Vtr (*-t) ^ x) dt 

0 


(4.46) 


In light of the discussion following Eq. (4.31), the error in Eq. (4.46) is expected to be 
small compared to the one-step solution, Eq. (4.31), for sufficiently small time steps. A 
truncated modal summation solution for this problem, similar to Eq. (3.13) for linear 
problems, is given by, 

. p 

T n = ^ 9n-kj z nj p « n (4.47) 

r*l 


or, in matrix form 

TUS n . k zi (4.48) 

Noting from Eq. (4.43) that the modes are now normalized with respect to K n . k and C n -k. 
the response given by the MAM is obtained by modifying Eq. (3.25) to yield 

i 

Tn = ^n-k Zn + 


-1 ~ 


Kn-k * ®n-k — ^n-k 
An-k 


Qn 


i-1 


(4.49) 


while that of the second-order FDM is obtained by modifying Eq. (3.35) as follows: 
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^ | J /s - *•> T j J 

Tn = <t>n-k Zn + K n -k * ^n-k ~ ^n-k Qn 


+ ' K n .k C n -k K n -k + ^n-k “ 4~ ^n-k On 

^ +• 


(4.50) 


The exact first-order time derivative used in Eq. (4.50) is obtained from Eq. (4.42) as 

• i*l 

For the purpose of easy implementation, Q„ has been approximated in this study as 


(4.51) 


<£‘a *u 


(4.52) 


. i-1 . 

The truncation of modes and the replacement of Qn by Ru arc the two numerical errors 
in the current form of the FDM. The former can be minimized by a proper selection of <l> 

**s 

based on the discussion in Sec. (3.3). The latter can be reduced by a frequent update of d> 
when necessary, that is, by maintaining a small value for k in Eq. (4.40). 

4.3 Modal Coordinates For a Piecewise Linear Time-Varying Load 

Recall that in Eq. (4.46), the nonlinear part of the generalized heat load vector 
is assumed constant during each time interval, and the integration is performed only for 
the explicit time-dependent portion. Accordingly, for a nonlinear problem when the 
linear part of the heat load vector, Rl. varies linearly with time Eq. (3.43) may be 
modified to yield, 

zh, - 2n.l, 4' V| “ * + 4W«( ( fc .1 - xir + (xft ) 


* Jfcr { (**» + QNli 4 ) 


(4.53) 
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The first two terms in Eq. (4.53) comprise the linear component of the modal coordinates 
while the third term represents the nonlinear contribution. 


4.4 Convergence Criterion and Distribution Error Norm 
The convergence criterion used to terminate the iteration process at any time t n is 
given by. 



£ E 


(4.54) 


where e is the specified tolerance. The distribution error norm of an approximation to the 
temperature vector is given by 


e = 


a a 

(T-T ) (T-T ) 


T 

TT 


(4.55) 


where T represents a full-system solution based on a well-refined mesh, and T 3 is an 
approximation based on the first p thermal modes. Note that this error norm can only be 
used to evaluate a method a posteriori, and is not intended to be used to predict the 
number of modes necessary for convergence of the modal solution. 


4.5 A Note on the Required Number of Modes 
The guidelines provided in Sec. 3.3 to estimate the number of modes required for 
linear problems cannot be applied directly to nonlinear problems. When the system 
matrices and hence the eigenmodes depend on the solution vector, which in turn depends 
on the number of modes included and the order of the method used, the issue becomes 
too complex for analysis. For linear problems, the improvement in the response obtained 
by the MAM as compared to the MDM is given directly by the pseudo steady-state term, 
CMAM, Eq. (3.51). For nonlinear problems, this linear superposition of the correction 
terms to obtain the response by the higher-order methods is not appropriate. It is 
therefore necessary to get the responses by the MDM, MAM, and FDM individually, 
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using their respective solution vectors to evaluate the nonlinear quantities. Then the 
correction made by the MAM, CMAM, is given by the difference between the entire 
right-hand sides of Eqs. (4.49) and (4.48), that is, by the difference between the responses 
of the MAM and the MDM. Similarly, CFDM is obtained from Eqs. (4.50) and (4.49) as 
the difference between the solutions of the FDM and the MAM. For the same reasons, it 
can be said that the superposition of the contribution of each successive higher mode to 
compute the total response of all the modes is permissible only for linear problems. 

4.6 A Note on the Computational Effort Involved 

The three modal methods presented in Sec. 4.2.2 for solving nonlinear thermal 
problems have been implemented in the COMET [63] on the CONVEX C220 high- 
performance computer at NASA Langley Research Center. The details of the 
implementation, which include the construction of the execution control file, the 
development of independent programs called processors, etc., following the rules and 
syntax of the COMET, can be found in Ref. [67]. The step-by-step procedure to be 
followed in the solution process using the modal methods is illustrated in Appendix A. 
The computational effort involved in the different steps for a nonlinear problem and the 
simplifications that occur for a linear problem are summarized below. 

The basic advantages of using the modal methods to solve thermal 
problems are: 

1. The need to solve much fewer equations than in the full system. 

2. The modal equations are uncoupled. 

3. The use of the convolution integral form of the solution for the modal coordinates. 

The skyline form of the symmetric conductance matrix, K, is used which requires less 
computer storage and also less computing time for all matrix operations involving K. 
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The main computationally intensive step for nonlinear problems is the solution of the 

EVP for intermediate basis updates, the frequency of which is problem-dependent. The 

product term K„!k On' 1 (Eqs. 4.49 and 4.50) involves the inverse of K. It is computed, 

however, without actually inverting the K matrix but by solving the linear set of 

equations K n -kr = Qn 1 whose solution, r, yields K n -kQn directly. Similarly, the term 

Kn-*k C Kn-k On 1 (Eq. 4.50) involves the solution of two sets of linear equations and a 

. i-1 

matrix-vector multiplication. Note also that in this study, Qn has been approximated as 

~ i-i 

r, in the EDM. The evaluation of Q^ 1 , (Eq. (4.41)) and the products d>^<X> Qn 
^ A 

(Eqs. 4.49 and 4.50) and S-L-C^Oi* 1 (Eq. 4.50) involve two matrix-vector products 

A 

each. Conventional implicit algorithms need to form and factor the left-hand side 
coefficient matrix and solve simultaneous algebraic equations at every time step for 
nonlinear problems. Although explicit algorithms avoid the soluriorf of simultaneous 
equations, they require much smaller time steps so that the solution remains stable. 
Modal methods, however, entail only the forward reduction and backward substitution 
stages in the solution of simultaneous equations, as the K matrix needs to be factored 
only with every eigenmode update and is held constant until the next update. The 
forward reduction and backward substitution processes could amount to significant 
computational effort though, when the solution is obtained in a large number of time 


steps. 

For the FDM, qA' 1 and Qn 1 can be treated as multiple loads and hence K n -k Qn 
and Knlk On 1 can be computed in parallel. Also for any method, the modal loads 


<lJ-kj qA 1 as well as the various matrix-vector products are good candidates for 
fine-grain parallel processing. 


For linear problems, the procedure is much simpler. Since all system matrices are 
assumed time-independent for the linear case, only one eigensolution and one 
factoriiation of K need to be performed. The guidelines presented in Sec. 3.3 can be 
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used to solve only for the desired number of eigenmodes. No iterations are required. 
Furthermore, if the applied load is constant or varies linearly in time at a constant rate, 
then no time-stepping is required either, as discussed already in Sec. 3.4, unlike nonlinear 
problems. Also, for the linear case, the response of the MDM and the MAM are obtained 
directly as by-products from the FDM. Nonlinear problems necessitate an independent 
run for each modal method, so that the response of that particular method is used to 
compute the temperature-dependent quantities which in turn affect the response. 



Chapter 5 

LINEAR TRANSIENT EXAMPLE PROBLEMS 

Literature survey indicates that the MDM is unsuccessful in the efficient solution 
of transient linear thermal problems, since excitement of higher frequencies in the wide 
spectrum response necessitates the inclusion of almost all modes for an accurate solution 
[54]. On the other hand, by representing approximations of the higher eigenmodes with 
few vectors in the basis, the Lanczos vectors are effective in the reduction process [56]. 
From the derivation of the higher-order modal methods presented in chapter 3, it is clear 
that these methods include an increasingly better approximation of the higher modes 
neglected by the MDM. If is therefore expected that the MAM and the FDM are more 
effective reduced-basis methods than the MDM for transient thermal problems and has 
been demonstrated for a linear one-dimensional problem [50]. 

This chapter introduces the application of the modal methods to linear thermal 
analysis before embarking on nonlinear analysis. In the following numerical study, the 
methods are evaluated based on the efficiency which is measured by the number of 
modes necessary to represent an accurate response. The accuracy of a solution is 
determined based on the distribution error norm, Eq. (4.55). Since the applied load varies 
linearly in time, the exact modal coordinates given in Eq. (3.40) were used to predict the 
response at any time in one step. 

5.1 Rod Subject to Convection at One End 

This simple linear transient example was chosen to introduce the modal methods 
to thermal problems and to demonstrate the improved efficiency of the higher-order 
methods compared to the MDM. The rod shown in Fig. 5.1 is initially at a uniform 
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Time, seconds 


A=a01 in 2 
L=Q1 in 
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k = 0.000214 BoVsec- in- °F 


Fig. 5. 1 Thermal finite element model of a rod subject to transient convective heating 
at one end (Linear example problem 1). 



53 


temperature of 0* F. Beginning at time zero, the left end is heated by convection while 
the right end is maintained at 0* F. A total of twenty elements was used to discretize the 
model. A forcing function that varies linearly with time was selected so that up to the 
second-order FDM could be applied. 

Temperature distribution at an early time t = 0.02 sec (when the gradient is still 
very steep) was obtained by the three modal methods (Figs. 5.2a - 5.2c). The response 
predicted by the MDM is shown in Fig. 5.2a. The full-system solution, denoted by the 
solid line in Fig. 5.2a, results from using all twenty modes. The acceptable error 
corresponds to a distribution that is fairly close to the exact solution. The solution 
approximated by using a reduced number of five modes underpredicts the peak 
temperature (417* F at left end) by 49%, and the distribution along the rod is also very 
oscillatory in nature. When the number of modes is increased to 10 and 15, the peak 
temperature is underpredicted by 23% and 11% respectively. Eighteen out of the 20 
modes are required by the MDM to yield a response with an acceptable error, Eq. (4.55), 
of 0.083. 

Figure 5.2b shows that the MAM overpredicts the peak temperature by 18% when 
using three modes and by 9% with four modes. However, with five modes the error falls 
within 0.08 and the oscillations in the distribution along the rod are mild. This faster 
convergence of the MAM compared to the MDM is due to the pseudo steady-state term 
which includes a first-order approximation of the higher modes that are neglected 
by the MDM. 

A higher rate of convergence is exhibited by the FDM which has the highest error 
with one mode but captures the distribution accurately with as few as three modes, 
Fig. 5.2c. It is evident that the derivative-related term in the FDM has further improved 
the approximation of the higher modes compared to the MAM. The rates of convergence 
of the three modal methods for this linear example with a C 1 continuous applied load are 
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x, inches 

a) Mode-displacement method (MDM). 

Fig. 5.2 Temperature distribution along a rod subject to transient convection at one end 
at time t * 0.02 sec (Linear example problem 1). 
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compared in Fig. 5.3. While the FDM requires just 15% of the modes (3 out of 20), 
which means a reduction of the problem size by 85%, the MAM and the MDM have 
achieved a reduction in problem size of 75% and 10%, respectively. 

5.2. Plate Subject to Uniform Surface Heating 

This problem was chosen to demonstrate the reliability of the measures presented 
in Sec. 3.3 to predict the number of modes required for linear problems. 

The plate shown in Fig. 5.4 is subject to a uniform specified heat load which 
increases linearly with time. The initial temperature of the plate is 0* F. All four edges of 
the boundary convect to the environment which is at a temperature T c = 0* F. 

The full plate was analyzed although this is a symmetric problem and only a 
quarter of the plate needs to be considered. The finite element mesh used for the modal 
analysis consists of 165 DOF with 140 uniform rectangular elements. This mesh is 
considered adequate for this problem since the full-system solution compares very well 
with that of a much finer mesh with 561 DOF. Since the load varies linearly with time, 
the modal solution at any time can be computed in one step. 

The guidelines listed in Sec. 3.3 are used to estimate a priori the number of modes 
that are required by the MDM, MAM, and the FDM at two different times t = 2 sec and 
9 sec, respectively. The norm of CMAM and CFDM were used to obtain the results 
presented herein. 

Estimate for the Required Number of Modes at Time t s 2 SCC 

The mode participation factor, Eq. (3.44), which indicates the component of the 
applied load that contributes to the response of the corresponding eigenvector, provides a 
measure of the significance of that vector in the total response. The participation factor 
of each mode is used to compute the error in the load vector, ei, Eq. (3.46), for increasing 
number of modes in the subset as shown in Fig. 5.5. Likewise, the errors, C2 and e 3» 


Number of inodes 


Convergence of the modal methods at time t * 0.02 sec (Linear example 

problem 1). 
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T(x,y,0) = 0.0°F 
= 0.02 Btu-in/in 2 -scc-°F 
k y = 0.01 Btu-in/in 2 -sec-°F 
p= 1.0 lb/in 3 
c= 1.0 Btu/lb-°F 


q, = 200t Btu/in 2 -scc 
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Fig. 5.4 Thermal finite element model (165 DOF) of a plate subject to uniform, 

transient heating over the surface and convection along the entire boundary 
(Linear example problem 2). 
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given by Eqs. (3.52) and (3.55), respectively, are computed for each additional mode and 
are shown in Figs. 5.6 and 5.7. The error converges in a step-like manner as seen in 
Fig. 5.5 indicating no change due to the addition of the second and third modes and then 
drops sharply with the inclusion of the fourth mode. This step-like convergence occurs 
because the second and third modes are asymmetric and therefore orthogonal to the 
uniform load and hence produce a negligible modal load. The step-like convergence over 
the entire spectrum is listed in Table 5.1. The wide range of modes 25 to 148 have 
negligible participation, but the 149th and the 161st modes cause a dramatic reduction in 
the error, thus confirming the importance of the higher modes in the MDM solution. 


Table 5. 1 A priori estimate of the required number of modes at time t = 2 sec 

(Linear example problem 2) 


Number 

of 

modes 

Error 

ei 

Number 

of 

modes 

Error 

C2 

Number 

of 

modes 

Error 

®3 

1 

MgEmm 

1 

0.598 E-02 

1 

0.383 E-03 

90 

0.443 E-01 

12-18 

0.138 E-03 

4-7 

0.452 E-04 

149 

0.261 E-01 

19-23 

0.126 E-03 

8 

0.578 E-05 

160 

0.246 E-01 

24 

0.927 E-04 



161 

0.576 E-02 

25-27 

0.284 E-04 




On the other hand, the error e2, Eq. (3.52), in the representation of the pseudo 
steady-state response, K* 1 R(t), also converges in a step-like manner. Fig. 5.6, but 
decreases rapidly due to the increasing magnitude of the higher eigenvalue appearing in 
the denominator. This also explains the faster rate at which the error e3, Eq. (3.55), 
approaches zero (Fig. 5.7) due to the square of the eigenvalue appearing in the 
denominator. The number of modes required by each of the methods is determined when 
the error approaches zero. From Figs. 5.5 - 5.7, it is estimated that up to 160 modes are 
required by the MDM whereas only 25 and 8 modes are required by the MAM and the 
FDM respectively. Referring to Table 5.1, these modes correspond to errors 
approximately two orders of magnitude less than the maximum error in each case. 
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Number of modes 

Fig. 5.6 Convergence of the error, e2, in the pseudo steady-state response used to 
° predict the number of modes required by the MAM at times t = 2.0 sec and 9.0 
sec (Linear example problem 2). 
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Fig 5 7 Convergence of the error, e 3 , in the derivative form of the pseudo steady-state 
response used to predict the number of modes required by the FDM at times 
t = 2.0 sec and 9.0 sec (Linear example problem 2). 
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Table 5.1 highlights the selected modes which cause the rapid convergence of the 
methods. 

Convergenc e of the Modal Methods at Time t = 2 sec 

The modal solutions obtained using the estimated number of modes at 
time t = 2 sec are shown in Figs. 5.8a - 5.8c. To judge the reliability of these estimates, in 
the case of the MAM for example, the distribution along y = 0 in. obtained by 10, 23 and 
25 modes are compared in Fig. 5.8b, since modes 13 through 23 are nearly orthogonal to 
the load. Similarly for the FDM the solutions from using eight and four modes are 
compared in Fig. 5.8c. Figures 5.8b and 5.8c demonstrate that 25 and 8 modes are indeed 
required by the MAM and the FDM respectively, and thus confirm that the corresponding 
estimates are accurate; whereas Fig. 5.8a shows that the MDM yields an acceptable 
solution even with about 90 modes, which is much fewer than the 160 modes estimated. 
This deviation from the estimate could possibly be due to the fact that the estimate is 
based on the error in the load vector and not the error in the solution vector itself. 

Estimate for the Required Number of Modes at Time t = 9 sec 

The rates of convergence of the errors, ei, t2 and e3, remain the same at this time. 
This is because the error e2, for instance, is essentially the error in the representation of 
K’ 1 by the subset of modes used. Since the eigenmodes are constant for this problem, 
this approximation does not change with time. The correction offered by the MAM, 
CMAM (Eq. (3.50)), varies with time since the load varies with time, as seen in Fig. 5.9. 
Although CMAM decreases as more modes are included at a given time, the maximum 
value of CMAM (which occurs when only one mode is included) increases with time as 
R (t) increases monotonically. For the given nature of the load, it is obvious from 
Fig. 5.9 that the MAM requires fewer modes to converge at time t = 9 sec than at time 
t = 2 sec. The rate of convergence of the correction of the FDM, CFDM (Eq. (3.53)), 
does not change with time (Fig. 5.10), since the applied load has a constant first 


x, inches 

a) Mode-displacement method (MDM). 

Temperature distribution at time t * 2.0 sec along y = 0.0 in. (Linear example 
problem 2). 
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Fig. 5.9 Convergence of the correction offered by the MAM, CMAM, at times t - 2.0 
sec and 9.0 sec (Linear example problem 2). 
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Number of modes 

Fig 5.10 Convergence of the correction offered by the FDM, CFDM, at times 
t = 2.0 sec and 9.0 sec (Linear example problem 2). 
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derivative with respect to time. Nevertheless, since the magnitude of CFDM is not 
negligible when compared to CMAM at time t = 9 sec, the FDM is expected to play an 
important role in the reduction process at this time also and hence converge with fewer 
modes than the MAM, albeit the difference is smaller at time t = 2 sec. The rate of 
convergence of the load vector, Eq. (3.45), is also expected to vary with time by realizing 
that the participation factor of the non-orthogonal modes is proportional to the load R(t). 
Although the corrections indicate that each method individually converges faster at this 
time than at time t = 2 sec, additional work may be necessary to use these corrections to 
decide the cutoff value for the number of modes required at different times. 

Convergence of the Modal Metho ds at Time t = 9 sec 

The modal responses were computed at time t = 9 sec and it is found that the 
MDM, MAM, and the FDM converge with fewer modes at this time, namely 30, 8 and 1 
respectively. The distributions shown in Fig. 5.11 or, alternatively, the distribution error 
norms compared in Fig. 5.12, confirm that all three methods converge with fewer modes 
at this time. For this problem with a monotonically time-varying load, the number of 
modes predicted based on the errors at time t * 2 sec is a conservative estimate, that is, 
guarantees convergence of the modal methods at all times but may not be necessary at all 
times, as the results at time t = 9 sec have shown. 

The results demonstrate the potential of the higher-order methods to effectively 
inprove the reduction achieved for transient thermal problems. Also the error estimates 
and the convergence of the correction terms show the potential to serve as useful tools in 
the prediction of the number of modes required at any given time. 
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Fig. 5.11 Temperature distribution at time t * 9.0 sec along y = 0.0 in. (Linear example 
problem 2). 



Chapter 6 

NONLINEAR TRANSIENT EXAMPLE PROBLEMS 

To study the feasibility of using the modal methods as a reduction technique for 
nonlinear transient thermal problems, two nonlinear examples were analyzed. The results 
of these analyses are presented in this chapter. The relative performance of the methods 
is evaluated and the conditions that bring out their best performances are highlighted. 
The multi-step approach described in Sec. 4.2.2 was used even though the load varies 
linearly in time, to accommodate for the temperature-dependent load vectors and 
eigenmodes. Also, the modal coordinates given by Eq. (4.53) were used to solve the 
numerical examples that follow. 

6.1 Rod Subject to Convection at One End 
with Temperature-Dependent Thermal Conductivity 

The problem statement for this simple example is shown in Fig. 6.1. It was 
primarily chosen to validate the new finite element algorithm for nonlinear problems 
which is presented in chapter 4. Two cases with different thermal conductivities are 

considered. 

6.1.1 Case l:k(T) = 0.0001 +0.5 E- 06 T 
Time for Eig en solution Update 

Some trial runs were made to determine the time when an eigensolution update is 
necessary to adequately represent the change in the nonlinear basis vectors with time. 
The frequency of these updates requires a compromise between computation time and 
accuracy. The FDM was used to march the solution from time t = 0 sec when the initial 
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Fie 6 1 Thermal finite element model of a rod subject to transient convective heating 

‘ at one end with temperature-dependent thermal conductivity (Nonlinear 

example problem 1). 
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eigenvalue problem was solved. An arbitrarily chosen time step, dt = 0.01 sec, 
(Eq. (4.36)) and a minimum of two modes (which is discussed subsequently) were used. 
The number of iterations needed to achieve convergence, within a tolerance of e = 0.01, 
increases with time as seen in Fig. 6.2, until the solution does not converge at time 
t = 0.08 sec, regardless of the number of modes in the subset. Even though the solution 
converges at time t = 0.07 sec, the number of iterations increases when more modes are 
included in the solution. A similar behavior was observed even when the tolerance was 
increased to e = 0.1, or, when the time step was decreased to dt = 0.005 sec. However, 
with an EVP update at time t = 0.05 sec, the number of iterations reduces to three from 
nine or ten at time t = 0.06 sec, and the increase in number of iterations at a later 
timet = 0.1 sec is less with more modes. This confirms the need to update the 
temperature-dependent eigensolution in time for nonlinear problems. However, this 
approach based on the iterations required is not meant to be used for an a priori estimate 
of the time for an EVP update. 

Convergenc e of Modal Methods 

The solution at time t = 0.1 sec in the early transient period when the gradient is 
significant, is used to evaluate the performance of the modal methods. A time step 
dt = 0.01 sec was used for this analysis with the initial eigensolution updated only once at 
time t = 0.05 sec, that is, m = 5 (see discussion following Eq. 4.41). The converged 
solution used to compute the distribution error norm of the reduced-basis solutions was 
obtained with a fine mesh of 101 nodes using a small time step dt = 0.001 sec. The 
implicit transient thermal analyzer in COMET which uses the Crank-Nicholson algorithm 
for time integration was used for this purpose. 

The MDM underpredicts the exact peak temperature of 373.7* F by 38%, 17% 
and 8.2% with 5 (25% of the number of degrees of freedom), 10 (50%) and 15 (75%) 
modes respectively. At least 16 modes are required by the MDM to yield a distribution 




Time, seconds 

Fig. 6.2 History of iterations required by the FDM to determine the time for an EVP 
update (Case 1 of nonlinear example problem 1). 
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fairly close to exact. On the other hand, the MAM converges with a much smaller subset 
of four modes while the FDM exhibits a superior performance by predicting the peak 
temperature within 1.36% using just two modes as shown in Fig. 6.3. This is made 
possible by the large ratio (10) of the magnitude of the derivative of the load to that of the 
load itself (which varies linearly with time) at this early time. The rates of convergence 
of the three methods are compared in Fig. 6.4. 

F.ffect of Time Step 

The effect of the time step size on the performance of the FDM was studied. The 
responses obtained by the FDM using two modes are compared for two different time 
steps dt = 0.01 and 0.05 sec in Fig. 6.5 which have an error norm of 0.038 and 0.08 
respectively. The smaller time step yields a distribution very close to the exact and helps 
to obtain a smooth temperature history when desired, as it computes the response at more 
instants of time. The advantages though, should justify the additional computational 
effort required to do the iterations at each of the intermediate times. Computations reveal 
that both time steps require the same total number of iterations, 17. The solution was 
further marched out in time to t = 0.2 sec. The solution at time t = 0.1 sec was used as the 
initial guess and the corresponding EVP was solved. The solution at time t = 0.2 sec 
obtained in one step (error = 0.09), is compared with that obtained in five steps using 
dt = 0.02 sec (error = 0.02) in Fig. 6.6. Again, the two analyses involve the same total 
number of iterations 7. These results are summarized in Table 6.1. By well-representing 
the transient variation of the nonlinear quantities, a smaller time step yields a better 
solution accuracy for a given computational effort. 

Minimum N umber of Modes 

It is seen in Fig. 6.4 that the error from the FDM is shown for two modes and 
more only. For the chosen time step (dt = 0.01 sec), tolerance (£ = 0.1) and EVP update 
(m = 5), it was determined by trial and error that a minimum of two modes were required 
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Fig 6 6 Effect of time step on the response of the FDM at time t = 0.2 sec with initial 
condition at time t » 0.1 sec and no EVP update (Case 1 of nonlinear example 

problem 1). 
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Table 6. 1 Effect of time step on the solution accuracy of the FDM 
(Case 1 of nonlinear example problem 1) 



Time = 0.1 sec 

Time = 0.2 sec 

Time step 

Number of 

Error 

Time step 

Number of 

Error 

dt 

iterations 

e 

dt 

iterations 

e 

0.01 

17 

0.038 

0.02 

7 

0.02 

0.05 

17 

0.08 

0.1 

7 

0.09 


by the FDM to obtain a solution (which is already converged), at t = 0.1 sec without any 
numerical instability or divergence at some intermediate time. A similar behavior was 
observed even with a fine mesh of 101 DOF. This minimum requirement for the number 
of modes did not arise for linear problems where the eigenmodes are independent of the 
solution and just one mode yields a solution even if highly erroneous. 

The time-marching parameters were varied to study their effect on the minimum 
number of modes required by the modal methods at t = 0. 1 sec. The results of this study 
are shown in Figs. 6.7a - 6.7c and are also summarized in Table 6.2. The MDM and the 
MAM yield a solution with one mode for all the parameters studied for this problem. The 
FDM also converges with one mode for the first set of parameters with only one EVP 
update and no iterations. Instead, when iterations are performed as with the second set of 
parameters, or, a piecewise linear approach (where the eigensolution is updated every 
time step, i.e., m = 1) is used with smaller time steps, the MAM converges with fewer 
modes (three instead of four) but the FDM requires more (two instead of one) for a 
minimum but the same number (two) to converge. 

6.1.2 Case 2: k (T) = 0.0001 + 0. 1 E - 04 T 

To study the effect of variation of the nonlinear parameter k(T), the slope of the 
conductivity curve was increased by a factor of 20. Although the heat load is identical in 
both cases, the diffusivity is very different. This is reflected in the temperature profiles at 
the same time t = 0.05 sec in Fig. 6.8. The temperature distribution in the first case has a 
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a) Mode-displacement method (MDM). 

Fig. 6.7 Effect of time-marching parameters on the number of modes required by the 

MDM, MAM, and the FDM at time t = 0.1 sec (Case 1 of nonlinear example 
problem 1). 
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relatively steep gradient compared to the second case, and the peak temperature is almost 
twice as well. 

Table 6.2 Effect of time-marching parameters on the number of modes required at time 
t - 0.1 sec (Case 1 of nonlinear example problem 1) 


Tii 

1 

ne -march 
parameter. 

ing 

5 

Minimum number 
of modes for a 
stable solution 

Modes required 
to converge within the 
desired accuracy 

dt 

m 

e 

MDM 

MAM 

FDM 

MDM 

MAM 

FDM 

0.01 

5 

1.0 

1 

1 

1 

17 

4 

2 

0.01 

5 

0.2 

1 

1 

2 

17 

4 

2 

0.01 

1 

1.0 

1 

1 

2 

17 

4 

2 

0.001 

1 

1.0 

1 

1 

2 

17 

3 

2 

0.0001 

1 

1.0 

1 

1 

2 

17 

3 

2 


Time for Eieensolution Update 

As in Case 1, the FDM was used to obtain the transient response starting with 
time time t = 0.0 as the initial condition. A time step dt = 0.001 sec was used which 
caused the solution to converge in only one iteration for a tolerance e = 1.0. A minimum 
of four modes were required in this case. As seen in Fig. 6.9, the error in the computed 
response progressively decreases from time t = 0 sec until it suddenly shoots up to 
prohibitive levels sometime after time t * 0.025 sec, even if more modes are included. 
With one EVP update at time t = 0.025 sec, the error at time t = 0.05 sec falls within an 
acceptable limit, confirming the need to update the nonlinear eigenmodes. 

Convergence of Modal Methods 

Using a time step dt = 0.001 sec and an EVP update at time t = 0.025 sec as 
determined above for the FDM, the solution at time t = 0.05 sec was obtained by the three 
methods. While the FDM required a minimum of four modes to converge, the MAM 
required seven. The MDM solution, though, became acceptable only with 16 modes. 
When the eigensolution was updated every 0.01 sec, the MAM also converged with only 





















































Time, seconds 

Fig. 6.9 Error history of the FDM to determine the time for an EVP update (Case 2 of 
nonlinear example problem 1). 
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four modes like the FDM. The rates of convergence of the three modal methods using 
dt = 0.001 sec and m = 10 are compared in Fig. 6. 10 for this case. 

F.ffect of Time Step 

For the given EVP update at time t = 0.025 sec, the unacceptable error resulting 
from using a time step dt = 0.025 sec is dramatically reduced by decreasing the time step 
to dt = 0.001 sec as seen in Fig. 6.11. Also, this entails little additional computational 
effort, as the solution requires only one iteration per time step. 

Minimum Number of Modes 

Results of a parametric study conducted similar to Case 1 are shown in 
Figs. 6.12a - 6.12c and are summarized in Table 6.3. For this case, the MDM again needs 
just one mode for a well-behaved solution with any of the parameters. With the first set 
of parameters, the MAM requires up to seven modes as a minimum and also for 
convergence. When the time step is reduced, this minimum decreases for the MAM, 
which then converges with four modes itself like the FDM. So for this case of higher 
nonlinearity, the reduction achieved by the MAM and the FDM is identical even at this 
early time t = 0.05 sec. On the contrary, the minimum requirement for the FDM 
increases when the time step is decreased, similar to the first case. 

In both cases of this example, the higher-order MAM and the FDM demonstrate 
superior convergence over the basic modal method, the MDM. Also the FDM is more 
effective than the MAM for a longer time in the first case with the less temperature- 
dependent conductivity. The difference in the convergence of the FDM relative to the 
MAM for the two cases is understood by studying the norm of the respective corrections 
obtained as discussed in Sec. 4.5. With two modes the proportion of the FDM 
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Fig. 6. 10 Convergence of the modal methods at time t = 0.05 sec (Case 2 of nonlinear 
example problem 1). 



Force-derivative method 
(Number of modes = 4) 


“ full-system 
■* dt = 0.001 sec 
~ dt = 0.025 sec 


Fig 611 Effect of time step on the response of the FDM at time t = 0.05 sec with initial 
g ‘ condition at time t = 0.0 and one EVP update at time t = 0.025 sec (Case 2 of 
nonlinear example problem 1). 
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a) Mode-displacement method (MDM). 

Fig 612 Effect of time-marching parameters on the number of modes required by 
the MDM, MAM, and the FDM at time t « 0.05 sec (Case 2 of nonlinear 
example problem 1). 






Table 6.3 Effect of time-marching parameters on the number of modes required at 
time t = 0.05 sec (Case 2 of nonlinear example problem 1) 
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Ti: 

1 

■ne -marching 
parameters 

Minimum number 
of modes for a 
stable solution 

Modes required 
to converge within the 
desired accuracy 

dt 

m 

£ 

MDM 

MAM 

FDM 

MDM 

MAM 

FDM 

0.001 

25 

1.0 

1 

7 

4 

16 

7 

4 

0.001 

10 

1.0 

1 

4 

4 

16 

4 

4 

0.001 

1 

10.0 

1 

4 

5 

16 

4 

5 

0.0001 

1 

10.0 

1 

2 

6 

16 

4 

6 

0.00005 

1 

10.0 

1 

2 

6 

16 

4 

6 


correction to that of the MAM looks significant at time t = 0.1 sec in Fig. 6.13a (Case 1), 
whereas in Fig. 6.13b (Case 2), with four modes it is almost negligible at 
time t = 0.05 sec. That is, CFDM contributes very little compared to CMAM, to the 
solution of Case 2. Figures 6.13a and 6.13b also illustrate that the relative importance of 
the higher-order terms changes with time. 

6.2 Lower Surface of Bay 3 of Shuttle Wing Segment 

The specific objective of this research is to compare the efficiency of the modal 
methods, in terms of the reduction achieved, in performing the nonlinear transient thermal 
analysis of a realistic structure, such as the Shuttle wing segment shown in Fig. 6.14. The 
thermal model shown in Fig. 6.15 represents a 58 in. segment of the lower surface of bay 
3 and consists of a 0. 1 19 in. thick aluminum sheet (to represent the structure) covered by 
a 1.36 in. thick layer of high-temperature reusable surface insulation (HRSI). The 0.16 
in. thick strain isolator pad (SIP) lies between the wing skin and the HRSI. The HRSI is 
bonded to the SIP and the SIP is bonded to the skin with room temperature vulcanized 
(RTV) adhesive. 


The lateral edges and the aluminum structure are assumed to be adiabatic. The 
outer wing surface is subject to heat loads representative of Shuttle reentry. The nature of 




















































Fig. 6.14 Geometry of Shuttle wing segment at Wing Station 240. 






100 


the heating history is significantly complex. As depicted in Figs. 6.16 and 6.17, in 
addition to the transient variation of the intensity, the spatial distribution of the load on 
the surface changes dramatically as well. This is a result of the air in the boundary layer 
over the wing undergoing transition from a laminar to a turbulent flow beginning at about 
700 sec into the heating history. Radiation loss from the surface is modeled using an 
equivalent nonlinear convection coefficient whose variation with temperature is shown in 
Fig. 6.18a. The thermal properties (specific heat and thermal conductivity) of the HRSI 
are highly temperature-dependent and their variation is shown in Fig. 6.18b. The various 
materials involved whose properties considerably differ between them, coupled with the 
complex imposed heating, result in a nonsymmetric temperature distribution. 

6.2.1 Simplified One-Dimensional Model 

In a preliminary analysis, the spatial variation in the applied heat load was ignored 
and uniform specified heating was assumed. This resulted in a simplified 
one-dimensional model across the thickness of the HRSI. Two-node conduction 
elements were used to model the HRSI, SIP, RTV and aluminum skin, and zero-length or 
point elements were used to represent the external heating and convection. The purpose 
of this analysis is to characterize the finite element discretization necessary across the 
HRSI thickness, so that the entire thermal response (not only the peak temperature) can 
be predicted with reasonable accuracy and minimal computational effort throughout the 
load history. Studies based on the simplified model are used to gain insight into the 
physics of this problem and to understand how it affects the performance of the modal 
methods. These preliminary results are expected to be useful in the subsequent finite 
element modeling of the actual two-dimensional problem. 

Mesh Convergence 

A series of meshes were generated starting with a relatively coarse mesh 
(Fig. 6.19) consisting of 16 nodes with ten uniform-length elements through the HRSI 
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Fig. 6. 16 Heating history for the lower surface of bay 3 of Shuttle wing segment (Also 
see Fig. 6.15). 


q*. Btu/in 2 - 










l.OE-4 


oC 

°. 8.0E-5 


8 

00 


£ 6.0E-5 


c 

o 


o 

s 

s 4.0E-5 


0 

1 


U 2.0E-5 I 


0.0E0 1 


0 


1000 


2000 


Temperature, °R 
a) Convection coefficient. 





Specific heat, Btu/lb-°R 


104 



Fig. 6.18 Concluded 
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b) HRSI thermal properties. 
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Fig. 6.19 Finite element mesh (16 DOF, uniform) for Shuttle one-dimensional thermal 
model. (Not to scale; all dimensions in inches). 
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thickness. Successively refined meshes include 26, 46, 80 and 166 DOF with 20, 40 50 
and 160 uniform elements, respectively, to represent the HRSI. Finally, based on the 
nature of the problem, the 16 DOF mesh with ten uniform elements was modified to have 
ten graduated elements with an arbitrary stretching factor of 2. This gives a large range 
of element sizes, with fine elements in the region near the heated end and very large 
elements in the interior. 

The full-system solutions were obtained using a relatively small time step 
dt = 0.5 sec with the eigensolution updated every time step. The well-refined 166 DOF 
system yields a very smooth solution as depicted by the solid lines in Figs. 6.20 and 6.21. 
The temperature at the HRSI surface rises very rapidly as the heating is applied and 
diffuses rather slowly through the HRSI thickness and SIP to the aluminum skin as 
shown in Fig. 6.20. After peak heating occurs at time t = 1 100 sec, the surface begins to 
cool while the temperatures of the interior of the HRSI and the aluminum skin continue to 
rise. The temperature distribution pattern in Fig. 6.21 shows steep gradients at early time 
due to the low conductivity of the HRSI. A good modal solution should capture this 
nonlinear behavior accurately with as few modes as possible. 

The 166 DOF solution is considered to be converged and is used to compute the 
distribution error norm of the other full-system solutions, to determine the adequacy of 
the meshes. As expected, the coarse 16 DOF systems both uniform and graded, have the 
highest error norm over the entire time period with a maximum of 0.008 and 0.017 
respectively at times t = 275 sec and t ■ 475 sec. The corresponding full-system solutions 
denoted by triangles and circles respectively (Fig. 6.21), are compared with the 
converged solution at discrete times when the errors are high. The distribution of the 
uniform mesh is not as smooth as those of the finer uniform meshes (not shown); 
nevertheless the nodal temperatures match those of the 166 DOF mesh at corresponding 
locations. The full-system solution of the graded mesh is far from smooth. The error at 
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Time, seconds 

Fig. 6.20 Temperature histories at different locations through a Shuttle rile and 
substructure obtained from a one-dimensional thermal model (Also see 
Fig. 6.19). 



Temperature, 
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y, inches 

Fig. 6.21 Comparison of full-system solutions of different meshes for Shuttle 
one-dimensional thermal model. 
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the location of the mid-HRSI thickness starts increasing around time t = 350 sec and 
reaches a maximum at time t = 475 sec (temperature distribution at intermediate times are 
not shown), when the peak temperature is almost reached but the temperature penetrates 
at a fast rate due to increase in conductivity, but the nodes are too far apart to capture this 
nonlinear behavior. Notwithstanding this error, the peak temperature as well as the 
gradient at the heated end are captured accurately at all times. The 46 DOF and 80 DOF 
solutions are very smooth and agree with the converged solution. 

Convergence of Modal Methods for the Various Meshes 

Having determined the limitations of the various meshes used and the accuracy of 
the respective full-system solutions, the modal solutions were obtained to identify the 
effect of the mesh on the eigensolution and hence on the efficiency of the modal methods. 
Although the transient variation of the load is curvilinear (Fig, 6.16), the thermal analyzer 
assumes linear variation between the specified input times so that the highest order modal 
method that can be used is the second-order FDM. For the purpose of a one-to-one 
comparison, a time step of 2 sec (determined by some trial and error) was used in all of 
the following analyses, (except where indicated otherwise) and the eigensolution was 
updated every 50 sec (that is, m = 25). 

With the crude 16 DOF uniform mesh, each of the three methods requires 12 
(75%) modes to converge, yielding a reduction of only 25%. While in the case of the 
graded mesh, the coarseness produces an inadequate solution accuracy in the interior of 
the HRS I, in the uniform mesh it decreases the efficiency of the MAM. The MAM 
converges with 20 modes on all other uniform meshes. This yields a steady increase in 
the percentage of reduction achieved with refinement beyond the 26 DOF mesh as seen in 
Fig. 6.22. It is noteworthy that with the 16 DOF coarse but graded mesh, the MAM 
shows a dramatic convergence with just 6 modes, 50% less modes than the uniform 
mesh, subject to the error bound of the full-system solution. This reduction is 
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comparable to that of the much finer 46 DOF uniform mesh. The increase in efficiency 
of the MDM is observed to be far less with the graded mesh which requires 1 1 modes 
(instead of 12 modes for the uniform mesh) and even with the finest 166 DOF mesh 
which still needs to retain 70 (42%) modes. Based on the requirement for convergence 
throughout the heating history, the FDM has not performed better than the MAM with 
any mesh (this issue is discussed later). Table 6.4 summarizes some of these results 
below. 


Table 6.4. Comparison of reduction achieved by the MDM and 
the MAM for Shuttle 1-D model 


Mesh description 

Required number of modes 

Reduction % 

MDM 

MAM, FDM 

MDM 

MAM, FDM 

16, uniform 

12 

12 

25.0 

25.0 

16, graded 

11 

6 

31.25 

62.5 

166, uniform 

70 

20 

58.0 

88.0 


Figure 6.23 compares the eigenvalues of the meshes at time t =* 350 sec and 
clearly shows how the mesh affects the distribution of the eigenvalues within the 
spectrum. The relative position of the fourth eigenvalue corresponding to the 16 DOF 
meshes, for instance, is given by 0.25 on the x-axis in Fig. 6.23. Mesh refinement 
enhances the rate of increase of the eigenvalues at the lower end of the spectrum, 
although the degree of enhancement decreases as the mesh becomes finer. A higher 
magnitude for the eigenvalues has been attained more efficiently, that is, by retaining 
fewer DOF in the mesh, by a suitable grading of the element sizes. This explains why the 
efficiency of the MAM with the 16 DOF graded mesh is comparable to the efficiency 
with the much finer 46 DOF mesh, and is much higher than that with the 16 DOF 
uniform mesh. A comparison of the mode shapes of the 16 DOF uniform and graded 
meshes shown in Figs. 6.24a and 6.24b shows how the steep gradients can be well- 
captured by fewer modes of the graded mesh. 
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Fig. 6.23 Comparison of eigenvalues at time t = 350 sec of the different meshes for 
Shuttle one-dimensional thermal model. 
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As in the previous examples, the MAM has achieved a remarkable improved 
reduction compared to the MDM for this one-dimensional model of the Shuttle example. 
For this example, although the heating history is transient, the FDM does not lead to any 
improvement over the MAM with any mesh. The number of modes required for 
convergence of the MDM and the MAM are based on the requirements for a stable as 
well as accurate solution over the entire time domain. 

The corrections CMAM and CFDM discussed in Sec. 4.5, were computed for the 
166 DOF mesh based on the initial eigensolution and the applied load at time t = 50 sec. 
As estimated from Fig. 6.25, the FDM and the MAM converge with four and eight modes 
respectively at time t = 50 sec and the corresponding distributions shown in Fig. 6.26 
confirm the accuracy of these estimates. The convergence of CMAM and CFDM at 
different times shown for the 166 DOF and 16 DOF graded meshes in Figs. 6.27a and 
6.27b were obtained using the corresponding full-system solutions at each time. As 
depicted qualitatively in these figures, the more-effective higher-order method and the 
modes required for convergence vary with time for this complex nonlinear problem with 
a transient load. That is, the correction offered by the FDM to the MAM becomes 
negligible rather early in time, whereas the correction offered by the MAM to the MDM 
remains significant even at a much later time t = 950 sec. This data was not used directly, 
however, to estimate the modes required for this problem. 

Effect of Transient Load and Nonli nearity on the Performance of the Modal Methods 

For this study a piecewise linear approach with dt = 0.5 sec was used. For the 
166 DOF mesh, the response of the FDM with four modes becomes unstable after 
time t = 200 sec. Both the MAM and the FDM converge with the same number of modes 
after time t = 250 sec. Similarly for the 16 DOF mesh, the MAM and the FDM converge 
with five modes up to time t = 300 sec. Recalling the implementation described in 
Sec. 4.6, the full form of the generalized load vector was used for the MAM but only the 



116 



Fig. 6.25 Convergence of the corrections offered by the MAM and the FDM, CMAM 
and CFDM respectively, at time t = 50 sec using the Shuttle one -dimensional 
thermal model with 166 DOF. 
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Fig. 6.26 Temperature distributions obtained by the MAM and the FDM at time 

t * 50 sec using the Shuttle one-dimensional thermal model with 166 DOF. 
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Number of modes 
a) 166 DOF uniform mesh. 

Fig. 6.27 Convergence of the corrections offered by the MAM and the FDM, CM AM 
and CFDM respectively, at different times for Shuttle one-dimensional 
thermal model. 
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approximate form of the derivative of the generalized load vector was used for the FDM. 
However, since the eigensolution was updated every time step in this study and no 
iterations were made, the nonlinear correction did not exist. Also for the 16 DOF graded 
mesh, the derivative of the nonlinear convection heat load vector (Fig. 6.28), computed 
based on the full-system solution, was included so that no approximation was actually 
made to the FDM. 

In the early stages, CFDM is significant and increases with time (Figs. 6.27a, 
6.27b) so that the FDM converges faster than the MAM. Later on, between 
times t = 250 sec and t = 300 sec, the trend reverses and CFDM decreases until the MAM 
becomes the most effective method. Figure 6.28 shows that the derivative of the applied 
heat load reaches a local peak at time t = 250 sec and then starts decreasing while the 
magnitude of the load itself is still increasing as seen in Fig. 6.29. The change in 
effectiveness of the FDM with time, which is experienced with both meshes, may be 
governed by the heating history for this problem. 

Although the MAM gives a converged solution with five modes for the 16 DOF 
mesh up to time t = 300 sec, the eigenvalues (Fig. 6.30), and hence the response, exhibit 
unrealistic oscillations in time. The difficulty faced by the MAM in predicting a stable 
solution after time t = 300 sec can be qualitatively explained with the help of Fig. 6.31. 
The normalized participation factor gives a measure of the contribution of each mode in 
representing the total load vector at different times. The spatial distribution of the load 
does not change with time for this 1-D model; however, the mode shapes vary with time 
due to the temperature-dependence of the thermal properties. Although the magnitude of 
the load is increasing with time, Fig. 6.31 shows that the participation of the modes 
steadily decreases from the initial time and reaches a minimum between times t = 300 sec 
and t = 400 sec when the inclusion of sufficient modes becomes necessary for a stable 
solution after this time. This pattern was observed with the finer meshes too, which gave 
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Fig. 6.28 Comparison of the histories of the first derivative of the specified and 
convective heat loads for the Shuttle one-dimensional thermal model. 
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Fig. 6.29 Comparison of the histories of the specified and convective heat loads for the 
Shuttle one-dimensional thermal model. 
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Fig. 6.31 Variation of the normalized mode-participation factors with time for different 
meshes of Shuttle one-dimensional thermal model. 
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a stable converged solution with 17 modes itself (instead of 20), provided the time step 
was decreased near time t = 300 sec. The MDM, however, is well-behaved even with 
much fewer modes than that required for convergence, although the solution is quite 
inaccurate as seen in Fig. 6.32. 

6.2.2 Two-Dimensional Model 

This example demonstrates the capability of the proposed higher-order modal 
method for solving complex situations with a much smaller set of equations than the 
MDM. Two-dimensional quadrilateral elements were used to model the HRSI, SIP and 
RTV while the aluminum structure, the heat load and the nonlinear surface convection 
were represented by line elements. Based on the experience with the 1-D model, the 
MAM was selected as the main focus in these analyses. For the 2-D model of the Shuttle 
problem (Fig. 6.15), the additional complexity involved is the spatial variation of the heat 
load. The magnitude of the load increases rapidly in a more or less spatially uniform 
manner from time t = 200 sec to time t « 350 sec (Fig. 6.17), and then the distribution 
becomes nonuniform. Also, around time t = 700 sec the pattern starts reversing, and by 
time t = 850 sec a step-like distribution pattern develops. Then the magnitude of the load 
increases at the maximum rate up to time t = 900 sec and after reaching the peak value at 
t = 1000 sec, the distribution becomes uniform again. A variation of the nonuniform load 
was also considered where the loads do not change abruptly as described earlier but in a 
much smoother manner as shown in Fig. 6.33. 

It is therefore prudent to arrive at the optimum mesh which can accurately 
represent the gradients, which in this case occur not only across the HRSI thickness but 
along the wing surface as well. For both cases of the spatially- varying load, judging from 
the 1-D analysis results, the discretization of the 16 DOF graded mesh was chosen across 
the HRSI thickness, i.e., along the y-axis in Fig. 6.15, as it offers a higher increase in 
efficiency for a given number of DOF. In an attempt to improve the solution accuracy, it 
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full-system; MAM, 6 modes 



Fig. 6.32 Comparison of temperature histories of the different modal solutions for the 
16 DOF graded mesh using a time step dt » 2 sec with the eigensoludon 
updated every 50 sec for Shuttle one-dimensional thermal model. 
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Fig. 6.33 Spatial distribution of heating on the lower surface of bay 3 at selected times 
(Smoothly-varying load). 



128 


was refined to include three more nodes at the mid-section where elements are rather 
large in the 1-D model. The discretization used along the x-axis is discussed separately 
for the two load cases. 

rase 1: Nonuniform Lo ad Distributio n ^ fh Tump Discontinuity 

The full-system solution at time t = 1000 sec obtained with a crude mesh of 187 
DOF shows unrealistic spatial oscillations as in Fig. 6.34, due to the improper 
representation of the applied load. This mesh was obtained using only two elements each 
to represent each different load span (Fig. 6.15). For this load case shown in Fig. 6.17, 
the discretization along the wing surface, i.e., along the x-axis in Fig. 6.15, must be done 
bearing in mind that the calculation of the nodal load vectors involves the Finite element 
approximation of lumping the element load equally at the two nodes. Consequently, a 
series of fine elements are necessary at the junction of the unequal input loads in order to 
accurately model the spatially-discontinuous load and thereby predict the temperature 
response reliably. The mesh thus generated after careful consideration consists of 578 
nodes as shown in Fig. 6.35. 

The temperature contours obtained from the full-system solution at 
time t = 1000 sec are shown in Fig. 6.36. It is evident that the temperature variation in 
the x-direction closely resembles the load distribution pattern causing a jump in the 
temperature of up to 150*R within a distance of 0.2 in. at locations where the load 
increases abruptly, and then the temperature remains constant (even if intermediate nodes 
are added) until the next jump location. To establish the suitability of the nodal locations 
along the x-axis for this problem, a mesh with 1683 nodes was employed along with a 
minor smoothing of the load changes over a thickness of about 1 in. The response was 
very similar to that of the 578 DOF mesh, that is, with similar jumps in temperature and 
temperature diffusion occurring only in the same width as the load smoothing in spite of 
the large array of very fine elements at the junctions. The adequacy of the discretization 
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Fig. 6.34 Temperature contours at time t = 1000 sec of Shuttle two-dimensional thermal 
model with 187 DOF and discontinuous load. 
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Fig. 6.35 Finite element mesh (578 DOF) for Shuttle two-dimensional thermal model 
with discontinuous load. 
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in the y-direction was verified by comparing the distribution along the edge x = 0.0 with 
the full-system solution of the 1-D model using the 166 DOF mesh. 

The convergence of the correction, CMAM, at time t = 350 sec (Fig. 6.37) is used 
to make a rough estimate of the number of modes required for convergence of the MAM. 
The reduced-system solution is considered to be converged when the error norm (based 
on the full solution) is acceptable throughout the solution time (Fig. 6.38). Based on this 
criterion, a subset of up to 400 modes (or 70%) are required by the MDM, whereas only 
200 modes (or 35%) are required by the MAM as seen in Fig. 6.39; thus confirming the 
improved predictability of the MAM for problems with extreme temperature gradients, 
even an abrupt increase in temperature caused by a discontinuous load. Here again, the 
FDM shows no better convergence than the MAM as expected from Fig. 6.37 which 
shows that CFDM is insignificant compared to CMAM. 

Case 2: Non uniform Lo ad Distribution with a Smooth Variation 

The sudden jump in the thermal load at discrete points in Fig. 6.17 was 
smoothened by incrementing the load in small step sizes over the entire load span and the 
resulting distributions at times t = 850 sec and 1000 sec are shown in Fig. 6.32. 
Accordingly, a vast number of elements were used to model these loads producing the 
mesh comprised of 986 DOF shown in Fig. 6.40. As expected, a smooth temperature 
distribution occurs along the wing surface as seen in Fig. 6.41, to obtain which 850 
modes (or 86%) are needed by the MDM (16% more than that required for Case 1). 
About the same increase in percentage of modes as that of the MDM is witnessed for the 
MAM which requires 550 modes (or 55%) for this case. 

The two cases discussed above, which encompass a wide range of nonsymmetric 
load distribution, firmly establish the superior performance of the MAM over the MDM 
although the degree of the increase in efficiency could vary, depending on the nature of 
the load distribution and the finite element discretization. 
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Fig. 6.37 Convergence of the correction offered by the MAM, CMAM, at time t = 350 
sec using the Shuttle two-dimensional thermal model with 578 DOF and 
discontinuous load. 



Error, e 


134 


MDM 



Fig. 6.38 Comparison of the error histories of the different modal solutions of the 
Shuttle two-dimensional thermal model with discontinuous load and 
578 DOR 
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Fig. 6.39 Convergence of the MDM and the MAM at time t = 1000 sec using the 
Shuttle two-dimensional thermal model with 578 DOF and discontinuous 
load. 
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Aluminum structure 



Fig. 6.40 Finite element mesh (986 DOF) for Shuttle two-dimensional thermal model 
with smoothly-varying load. (Also see Fig. 6.33). 
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Fig. 6.41 Temperature contours at time t » 1000 sec of Shuttle two-dimensional thermal 
model with 986 DOF and smoothly-varying load. 




Chapter 7 

CONCLUSIONS AND RECOMMENDATIONS 

The force-derivative method (FDM), which collectively represents a class of 
higher-order modal methods that improve the fundamental modal method, the MDM, 
with various-order correction terms, has been presented for nonlinear transient thermal 
analysis. The additional terms in the FDM include the forcing function and its 
derivatives with respect to time. A new algorithm incorporating the modal methods and a 
fixed-point iteration scheme has been implemented in an existing advanced finite element 
code, COMET, and validated with the aid of numerical examples. The solution is 
advanced in time with the nonlinear system matrices and load vectors being re-evaluated 
during the iterations at each time step, while the eigensolution is updated periodically to 
account for the change in the nonlinear basis vectors. For nonlinear problems, the first- 
order correction that the MAM offers to the MDM is fully realized by forming a 
generalized load vector, which in addition to the applied load vectors includes a 
corrective vector to account for the change in the nonlinear eigensolution between 
updates. The similar implementation of the second-order FDM in its entirety, requires 
the exact derivative of the generalized load vector which entails additional computational 
effort. Hence, in general, the highest-order FDM that may be employed, could be 
decided based on the order of the explicit time-dependence of the forcing function, as for 
linear problems. 

In general, the results demonstrate the potential of the higher-order methods to 
effectively improve the reduction achieved for transient thermal problems. In this study, 
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the efficiency of the methods is measured only by the reduction achieved, in terms of the 
number of modes required, to obtain the response within the desired accuracy. 

The first nonlinear example is a rod with a transient load and temperature- 
dependent thermal conductivity. The second nonlinear example involves a realistic 
structure, the lower surface of bay 3 of the Shuttle wing segment with nonlinear thermal 
properties. The complicated heat loads include specified surface heating which varies not 
only with time but in space as well, and nonlinear surface convection. 

A number of factors that affect the performance of the modal methods have been 
identified. The need for a multi-step approach in time with periodic updates of the 
eigensolution has been established for nonlinear problems. Reduction of the time step 
size and eigensolution updates have either increased the solution accuracy or improved 
the convergence of the methods to some extent. The magnitude of the derivative of the 
transient load compared to that of the load itself, and the degree of temperature- 
dependence of the thermal conductivity are seen to affect the relative effectiveness of the 
FDM compared to the MAM. 

The preliminary one-dimensional analysis of the Shuttle problem using different 
meshes clearly indicates that for a given number of degrees of freedom a suitably-graded 
mesh based on the expected response can upgrade the eigenmodes and thereby further 
enhance the faster convergence of the MAM over the MDM. For the two-dimensional 
problem with a discontinuous load distribution, the mesh is contrived in light of the above 
conclusion. Results confirm that the correction term of the MAM, which involves the 
load itself, is very effective in representing the neglected higher modes thus enabling the 
MAM to achieve a remarkable reduction of 65%, over twice that of the MDM, for this 
example with complex loading conditions. 

The above examples have demonstrated that the FDM (the MAM in particular) is 
indeed a feasible, effective reduction method for nonlinear transient thermal problems. 
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To solve nonlinear problems more efficiently, future research must include easy but 
effective means to approximate the changing eigensolution. When coupled with a 
transient adaptive meshing scheme, the method shows potential to yield a reliable and 
efficient solution for problems with severe gradients. 

In this study, it has been observed that the response of the MDM, though highly 
inaccurate, is well-behaved even with much fewer modes than that required for 
convergence. However, a minimum number of modes are required by the higher-order 
methods to yield a stable solution without unrealistic oscillations in time, regardless of 
how small the time step may be. This behavior of the MAM and the FDM exhibited in 
this study needs further investigation. 

Error estimates based on the approximation of the pseudo steady-state response 
and its time-derivative by a subset of modes have been identified for linear transient 
analysis. Results of a linear problem with a spatially-uniform load but a linearly-varying 
transient load, although not conclusive, show how these error estimates can reliably 
predict the number of modes required by the MAM and the FDM throughout the time 
domain. The results also indicate the potential usefulness of the convergence of the 
correction terms of the MAM and the FDM in determining the effective method and the 
modes required, both of which may vary with time. Further study is required for a 
thorough interpretation of these correction terms, to decide how they may be obtained for 
nonlinear problems and for establishing convergence criteria to use them as a priori 
estimates. 

The parallelization of the FDM is another subject identified for future study. The 
uncoupled nature of the modal equations and the numerous matrix-vector products 
involved in the modal solutions indicate that the computational efficiency can be 
improved by using parallel processing techniques. 
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Appendix 

COMPUTATIONAL PROCEDURES USED IN THE IMPLEMENTATION 
OF THE MDM, MAM AND THE FDM 

The sequence of the steps to be followed in the solution of transient thermal 
problems using the modal methods is presented below. 

1 . Input the number of modes to be included in the solution. 

Input the number of time steps before every EVP update, m. 

At initial time t = 0: 

Evaluate the nonlinear system matrices, Kq and Co- 

Evaluate the linear and nonlinear applied load vectors, Rlq and Rnlq, and their time 
derivatives, Rlq and RnLq- 

2. Apply the boundary conditions. 

A A 

3. Solve the initial EVP, Eq. (4.25), to get Oo and Ao. 

Factor Ko. 

Set the indicator for EVP update, k, to zero. See Eq. (4.40). 

4. For n = 1 , 2 . . . total number of timesteps; loop through step 27. 

Increment k; k = k+1. 

5 . Evaluate the linear applied load and its time derivative at current time, RLj, and Rl„- 
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6. Initialize the modal coordinates, z n -i, Eq. (4.45). 
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7 . Compute the linear part of the modal coordinates; sum of the first two terms on the 
right-hand side of Eq. (4.53). 

8. For i = 1, 2 . . . until the iterations converge, loop through step 24. 

9. Compute the nonlinear corrective load vector, Qnl^ 1 . Eq. (4.41), and add to the 
nonlinear applied load vector, Rnl^ 1 to form the nonlinear part of the generalized 
load vector. 

10. Obtain the modal coordinates, z n , by adding the nonlinear part (third term on the 
right-hand side of Eq. (4.53)) to the linear part computed in step 7. 

11. Get the approximate response, Tn, by the MDM, Eq. (4.48). If the MDM is the 
chosen method, go to step 21. 

12. Otherwise, assemble the generalized load vector Q^ 1 , Eq. (4.42), by adding the 
linear component from step 5 to the nonlinear component from step 9. 

1 3. Solve the set of linear equations K n -k r = Q^ 1 to obtain K n lk Qh \ 

14. Compute <DQn'\ 

A 

15. Get the response of the MAM, Eq. (4.49). If the MAM is the desired method, 
go to step 21. 
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16. 


. i-l . i-1 

Otherwise, using the approximate form of Q n , Eq. (4.52), solve K n -k r = Q n to 
yield K n ! k Qi \ 


17. Compute C n -k K n ! k OJ, \ 


18. Solve Kn-k r— Gn-k Kn-k Qn to yield Kn- k Cn-k Kn- k On • 


^ 1 ~ T .i-l 

19. Compute G> 0 Qn • 

A 

20. Get the response of the second-order FDM, Eq. (4.50). 

21. Use the new temperature vector to update the nonlinear system matrices and 
load vector for the next iteration or time step. 


22. Apply the boundary conditions. 


23. If the iterations have converged, Eq. (4.54), go to step 25. 


24. Otherwise, go to the next iteration, step 8. 

25. If this is the last time step, STOP. 

26. Otherwise, check if it is time to update the EVP; if it is not, i.e., k * m, go to the next 

/N /S 

time, step 4. If it is, then obtain the new set of eigenmodes, <D n and A n . 

Factor K„. 

Reset k = 0. 

Go to step 4. 
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